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Prediction  of  waves  and  currents  in  the  nearshore  region  in  the  presence  of  proposed 
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and  a  circulation  model.  The  wave  model  uses  a  parabolic  approximation  to  a  combined 
refraction-diffraction  mild  slope  equation  which  includes  currents.  The  circulation  model 
uses  an  alternating  direction  implicit  method  solution  to  the  depth  averaged  equations  of 
momentum  using  the  radiation  stress  terms  obtained  through  solution  of  the  wave  model  as 
the  forcing  terms.  Each  of  the  terms  in  the  depth  averaged  equations  of  momentum  involve 
assumptions  that  are  discussed  in  the  report.  The  model  is  constructed  in  modular  form  so 
that  improvements  can  easily  be  incorporated  into  the  model. 

The  model  was  employed  for  several  test  cases  such  as  rip  channels,  bar  gaps,  and 
both  shore  perpendicular  and  shore  parallel  breakwaters.  Good  agreement  was  obtained 
when  the  model  was  used  to  simulate  a  laboratory  experiment  of  Michael  R.  Gourlay 
("Wave  set-up  and  wave  generated  currents  in  the  lee  of  a  breakwater  or  headland," 
Proc.  14th  Coastal  Eng.  Conf.,  Copenhagen,  1974,  pp.  1976-1995). 
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Suggestions  are  made  for  improving  the  model  and  for  including  sediment  transport  in 
the  model. 


CHAPTER  1 
INTRODUCTION 


As  waves  propagate  they  are  transformed  by  variations  in  currents  and  depths,  along 
with  interaction  with  obstructions.  The  major  processes  of  wave  transformation  include 
shoaling,  reflection,  refraction,  diffraction  and  dissipation.  Shoaling  is  the  change  of  wave 
energy,  propagation  speed  and  wave  height  due  to  changes  in  the  water  depth  and/or  current 
field.  Refraction  is  the  change  in  the  direction  of  wave  propagation  due  to  the  spatial 
gradient  of  the  bottom  contours  and/or  current  field.  Diffraction  is  the  phenomenon  of 
wave  energy  being  transmitted  laterally  along  a  wave  crest.  Dissipation  is  the  process  of 
energy  loss  due  to  a  multitude  of  possible  causes  such  as  bottom  percolation,  damping 
due  to  viscous  effects,  turbulence,  and  wave  breaking.  As  waves  approach  a  shoreline  and 
break,  momentum  exchanges  can  produce  currents  and  mean  water  level  variations.  These 
currents  can  be  in  the  form  of  longshore  currents,  rip  currents,  or  circulation  cells  in  the  lee 
of  obstructions.  They  in  turn  modulate  the  incident  wave  field.  The  process  can  be  quite 
complicated  depending  upon  the  complexity  of  the  bathymetry,  the  ambient  current  field, 
and  the  presence  of  obstructions. 

Coastal  engineers  and  planners  have  many  reasons  for  wanting  to  obtain  accurate  pre- 
dictions of  waves  and  currents  and  water  elevations  in  the  coastal  zone.  For  instance,  pre- 
dictions of  extreme  conditions  are  needed  for  site  planning  and  erosion  mitigation  studies. 
Harbor  authorities  need  wave  predictions  for  structure  design  and  for  operational  decisions. 
This  dissertation  develops  a  numerical  model  for  obtaining  predictions  of  waves,  currents 
and  mean  surface  elevations  in  a  given  domain  with  the  offshore  wave  conditions  and  the 
bathymetry  as  the  input.  The  model  uses  a  linear  wave  equation  and  depth  integrated 
equations  of  momentum  and  mass  conservation.    Comparison  with  experimental  results 
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indicates  that  the  model  produces  useful  and  valid  predictions  of  the  waves  and  currents  in 
the  nearshore  region  in  the  presence  of  obstructions  such  as  groins  and  breakwaters.  The 
predictions  obtained  through  use  of  the  computer  model  herein  reported  are  not  intended 
to  be  the  sole  resource  for  the  engineer.  A  good  engineer  or  planner  will  seek  answers  from 
multiple  sources  and  cross  check  them,  recognizing  the  limitations  of  each  method  used  to 
obtain  the  answer. 

Predictions  of  waves  and  currents  can  be  obtained  for  different  cases  analytically, 
through  experiments  or  by  numerical  solutions.  Analytical  solutions  are  obtained  only 
for  simple  cases  involving  idealized  beaches  using  simplified  assumptions  which  are  not  al- 
ways satisfied  in  nature.  Longuet-Higgins  and  Stewart  (1963)  give  an  analytic  solution  for 
wave  set-up  for  steady  state  conditions  on  a  plane  beach  with  straight  and  parallel  con- 
tours, assuming  linear  waves  approaching  the  beach  normally,  and  spilling  breakers  where 
the  wave  height  shoreward  of  the  breaker  line  is  a  constant  times  the  total  depth  of  the 
water  column.  Longuet-Higgins  (1970a)  gives  a  cross  shore  profile  of  the  longshore  current 
for  a  plane  beach  with  straight  and  parallel  contours,  for  linear  waves  approaching  the  beach 
at  an  angle,  spilling  breakers,  and  linearized  bottom  friction,  while  neglecting  the  influence 
wave  set-up  has  upon  the  total  depth  of  the  water  column.  In  a  companion  paper,  Longuet- 
Higgins  (1970b)  extended  the  solution  to  include  lateral  diffusion  of  the  current  momentum. 
There  are  other  analytical  solutions  for  waves  and  currents  but  all  are  limited  to  very  simple 
cases. 

The  literature  is  replete  with  laboratory  investigations  of  waves  and  currents.  Physical 
models  have  advantages  and  disadvantages.  The  advantages  of  experimentation  are  that 
results  can  be  obtained  for  complex  situations  where  there  is  uncertainty  about  the  exact 
form  of  the  governing  equations  or  about  the  precise  form  of  any  of  the  component  terms 
of  the  equation.  Any  reasonable  amount  of  detail  can  be  included.  Visualization  of  the 
physical  processes  provides  an  opportunity  to  identify  problems  which  otherwise  might  not 
be  evident. 
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However  there  are  limitations  to  the  use  of  laboratory  experiments.  There  is  the  problem 
of  scale  effects  which  means  that  all  scaling  requirements  cannot  be  satisfied  simultaneously; 
thus,  all  features  may  not  be  interpreted  properly  to  the  prototype.  For  example,  the 
time  scale  dictated  by  Froude  modelling  laws  in  which  gravitational  forces  dominate  is 
different  from  that  obtained  from  Reynolds  modelling  laws  in  which  viscous  forces  are  most 
important.  Another  disadvantage  of  physical  models  is  the  cost  of  labor  to  construct  and 
run  large  scale  experiments.  Also,  because  of  the  cost  of  construction,  it  is  not  always  very 
easy  to  change  details  in  the  model.  Thus,  it  is  not  always  possible  to  isolate  one  parameter 
for  study.  For  example,  an  investigation  of  the  influence  of  the  size  of  the  gap  between 
segmented  breakwaters  might  require  many  reconstructions  of  the  rubble  mound  structure. 
The  large  size  of  a  model  means  that  it  is  not  easily  stored  once  the  initial  experiment 
is  completed.  The  options  are  to  either  dismantle  the  model  upon  completion  of  testing 
or  to  incur  some  maintenance  costs.  The  significance  of  this  is  that  the  model  cannot 
always  be  easily  retested  if,  for  example,  additional  input  information  becomes  available  or 
man  or  nature  makes  a  substantial  alteration  of  the  prototype  after  the  model  has  been 
dismantled.  Also  the  results  of  a  model  for  one  prototype  are  not  always  transferable  to 
other  prototype  situations.  Often  the  available  space  for  the  physical  model  will  require 
the  use  of  a  distorted  model  in  which  the  vertical  scale  is  different  from  the  horizontal 
scale.  This  presents  the  additional  problem  that  for  a  distorted  model,  when  modelling 
a  situation  where  both  refraction  and  diffraction  are  equally  important,  the  laws  of  scale 
modelling  dictate  different  time  scales  for  refraction  and  diffraction.  For  refraction  the 
time  ratio  is  Tr  =  \/Lvr  where  the  subscript  r  is  for  the  ratio  of  the  value  in  the  model 
to  the  value  in  the  prototype  and  T  and  Ly  are  time  and  vertical  length  respectively.  For 
diffraction  Tr  =  i"t  in  shallow  water  and  TT  =  \JIjr,  in  deep  water,  where  Lh,  represents 
the  horizontal  length  scale. 

Numerical  models  do  not  have  many  of  the  disadvantages  found  with  physical  models. 
Any  mechanism  for  which  the  governing  equations  are  known  can  be  modeled  numerically. 
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There  are  no  limits  in  terms  of  scale  effects.  Numerical  models,  once  developed,  can  be 
adapted  easily  to  many  different  prototype  situations.  There  is  no  problem  in  the  storage  or 
duplication  of  the  model.  Today  computer  facilities  are  far  more  available  than  laboratory 
facilities.  Computer  models  easily  allow  for  sensitivity  testing  of  the  importance  of  an 
individual  variable  and  any  specific  feature  of  the  model  can  be  monitored. 

Numerical  models  also  have  their  disadvantages.  In  fact  some  of  their  advantages 
are  also  their  disadvantages.  For  example,  being  able  to  monitor  any  feature  also  means 
that  the  computer  will  deliver  only  the  information  that  is  explicitly  specified.  Thus  an 
investigation  of  one  aspect  will  not  indicate  problems  with  other  aspects.  Though  any 
mechanism  for  which  the  governing  equation  is  known  can  be  modeled,  the  model  is  only  as 
good  as  the  assumptions  involved  in  the  equations.  For  example,  the  forces  exerted  on  the 
water  column  by  the  complex  interaction  between  the  fluid  motion  and  the  bottom  are  not 
totally  understood  and  are  often  greatly  simplified.  Another  example  is  the  use  of  depth 
integrated  two-dimensional  models  of  fluid  flows  when  the  reality  is  three-dimensional.  Both 
of  the  above  are  incorporated  in  the  present  study.  Knowing  the  correct  equation  does  not 
always  make  for  successful  computer  modelling  since  the  complexity  of  the  equations  may 
make  for  difficulties  (for  example  solving  for  the  complete  three-dimensional  flow  field  using 
the  Navier-Stokes  equations).  Computer  time  and  cost  for  large  programs  can  often  be 
significant. 

1.1     Literature 

This  section  is  a  brief  history  and  literature  review  of  the  attempts  of  the  coastal 
engineering  profession  to  obtain  predictions  of  waves  and  currents  in  the  nearshore  region 
given  the  offshore  wave  conditions. 

The  extension  of  Snell's  law  which  governs  optical  wave  refraction  to  the  analogous 
problem  of  water  waves  over  straight  and  parallel  contours  has  been  known  for  some  time. 
Johnson,  O'Brien,  and  Isaacs  (1948)  reported  graphical  methods  for  constructing  wave 
refraction  diagrams.  Two  methods  were  reported,  one  producing  wave  fronts  and  the  second 
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producing  wave  rays.  Pierson  (1951)  showed  that  the  method  using  wave  fronts  introduced 
more  error  than  the  wave  ray  method.  These  methods  proceeded  upon  the  assumption  that 
the  contours  are  locally  straight.  The  wave  rays  obtained  were  important  in  that  the  wave 
height  was  determined  to  vary  inversely  as  the  square  root  of  the  spacing  between  the  wave 
rays.  Using  these  ray  construction  techniques  O'Brien  (1950)  was  able  to  demonstrate  that 
the  April  1930  destruction  of  the  Long  Beach  Harbor  breakwater  on  a  seemingly  calm  day 
most  probably  resulted  from  the  focusing  of  long  wave  energy  from  a  South  Pacific  storm. 
Munk  and  Arthur  (1952)  showed  a  method  to  obtain  the  wave  height  directly  from  the 
curvature  of  the  ray  and  the  rate  of  change  of  the  depth  contours.  The  wave  ray  tracing 
techniques  did  not  include  the  effects  of  currents  and  had  the  disadvantage  (from  a  numerical 
modeling  point  of  view)  that  wave  heights  and  directions  were  obtained  at  positions  along 
a  ray  and  not  at  pre-specified  grid  positions.  Noda  et  al.  (1974)  solved  this  problem  by 
devising  a  numerical  scheme  that  solved  a  wave  energy  equation  and  a  statement  for  the 
irrotationality  of  the  wave  number  to  obtain  wave  height  and  direction  at  grid  points.  This 
scheme  also  included  the  effects  of  current  refraction  but  neglected  diffraction. 

Numerical  models  of  wave  induced  currents  were  developed  by  Noda  et  al.  (1974)  in 
which  time  dependent  and  advective  acceleration  terms  were  ignored.  Birkemeier  and  Dal- 
rymple  (1976)  wrote  a  time-marching  explicit  model  which  neglected  lateral  mixing  and 
advective  acceleration.  Ebersole  and  Dalrymple  (1979)  added  lateral  mixing  and  the  non- 
linear advective  acceleration  terms.  Vemulakonda  (1984)  wrote  an  implicit  model  which 
included  lateral  mixing  and  advective  acceleration  terms.  All  of  the  above  used  the  numer- 
ical scheme  of  Noda  et  al.  (1974)  to  obtain  wave  heights  and  wave  angles.  Each  of  these 
models  neglected  diffraction. 

Diffraction  of  water  waves  represents  some  very  difficult  mathematical  problems.  His- 
torically most  solutions  were  based  upon  the  assumption  of  a  flat  bottom,  since  to  quote 
Newman  (1972), 

The  restriction  to  constant  depth  h  is  important Problems  involving  wave 

propagation  in  a  fluid  of  variable  depth  are  of  great  interest,  and  the  resulting 


diffraction  effects  are  among  the  most  important  of  all  ...  but  these  problems 
are  less  tractable  than  those  involving  a  fluid  of  constant  depth,  (page  2) 

Penney  and  Price  (1944)  showed  that  Sommerfeld's  solution  of  the  optical  diffraction  prob- 
lem as  described  by  Bateman  (1944)  is  also  a  solution  of  the  constant  depth  water  wave 
diffraction  problem.  Penney  and  Price  (1952)  solved  the  constant  depth  diffraction  problem 
for  breakwaters  assuming  a  semi-infinitely  long,  infinitesimally  thin  barrier.  The  solution 
is  in  terms  of  Fresnel  integrals  which  need  to  be  evaluated  numerically.  The  Shore  Protec- 
tion Manual  (SPM,1977)  produced  many  graphical  solutions  based  upon  Penney  and  Price 
(1952)  that  give  the  wave  height  and  direction  in  the  lee  of  breakwaters. 

The  problem  of  combined  refraction  and  diffraction  according  to  the  SPM  could  be 
handled  by  assuming  a  constant  depth  for  a  few  wavelengths  so  as  to  obtain  the  diffraction 
effects  and  then  continue  with  the  refraction  effects.  Perlin  and  Dean  (1983)  used  essentially 
the  same  procedure  in  a  numerical  model  which  calculated  the  longshore  thrust  in  terms  of 
the  breaking  wave  height  and  angle. 

Liu  and  Mei  (1975)  solved  for  the  wave  field  around  shore  parallel  and  shore  perpendicu- 
lar breakwaters  by  using  a  curvilinear  coordinate  system  where  one  coordinate  corresponded 
to  the  wave  ray  and  the  other  to  a  curve  of  constant  phase.  Their  solution  was  for  "an  am- 
plitude modulation  factor  D(£,$)  which  accounts  for  diffraction"  (page  72),  where  £  and  f 
are  the  two  coordinates.  An  equation  for  Z7(£,f)  was  obtained  and  a  boundary  value  prob- 
lem was  established  by  matching  conditions  in  the  incident  zone  and  the  shadow  zones.  For 
constant  depth  the  solution  for  Z?(£,f)  is  in  terms  of  Fresnel  integrals.  Upon  obtaining  the 
wave  field  the  radiations  stresses  were  obtained  and  used  to  force  a  circulation  model  that 
solved  for  the  stream  function  and  mean  water  surface  level  while  ignoring  the  non-linear 
advective  acceleration  terms  and  the  lateral  diffusion  of  current  momentum.  Although  Liu 
and  Mei  produced  results  for  the  wave  field,  currents,  and  set-up  resulting  from  the  diffrac- 
tion of  waves  by  breakwaters,  their  study  is  not  a  wave-current  interaction  model  since 
there  are  no  current  terms  in  the  wave  equations  and  no  feedback  of  current  terms  from  the 
circulation  model. 
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Berkhoff  (1972)  derived  an  equation  governing  combined  refraction  and  diffraction  with 
the  assumption  of  a  mild  slope.  Radder  (1979)  applied  the  splitting  method  of  Tappert 
(1977)  to  obtain  a  parabolic  equation  governing  the  propagation  of  the  forward  propagat- 
ing component  of  the  wave  field.  Booij  (1981)  extended  the  work  of  Berkhoff  to  include 
the  effects  of  currents  and  also  obtained  a  parabolic  approximation.  Kirby  (1984)  made 
corrections  to  the  equation  of  Booij . 

Parabolic  wave  equations  have  recently  been  used  in  wave-induced  circulation  models. 
Liu  and  Mei  (1975)  solved  for  the  wave  field  by  matching  conditions  between  the  shadow 
zones  and  the  illuminated  zone  in  the  lee  of  a  structure.  The  wave  information  was  then 
used  to  force  a  circulation  model  which  formulated  the  governing  depth  integrated  time 
averaged  momentum  equations  in  terms  of  a  stream  function.  However  this  model  did  not 
have  current  feedback  in  the  wave  equation.  Yan  (1987)  and  the  present  study  both  include 
a  current  feedback  for  the  parabolic  wave  equation  and  solve  directly  for  the  mean  currents. 

Prior  to  the  use  of  parabolic  wave  equations  to  obtain  the  wave  forcing  in  circulation 
models,  refraction  wave  models  using  a  numerical  scheme  derived  by  Noda  et  al.  (1974) 
were  used.  Birkemeier  and  Dairy mple  (1976)  produced  an  explicit  numerical  solution  of 
the  depth  averaged  equations  of  momentum  and  continuity  using  the  wave  model  of  Noda 
et  al.  (1974).  Their  model  did  not  include  the  nonlinear  advective  acceleration  terms 
or  any  mixing.  Ebersole  and  Dalrymple  (1980)  included  the  advective  acceleration  terms 
and  lateral  mixing  in  an  explicit  model  which  also  used  the  refraction  scheme  of  Noda 
et  al.  Vemulakonda  (1984)  Produced  an  implicit  model  which  also  was  based  upon  the 
wave  refraction  scheme  of  Noda  et  al.  Yan  (1987)  used  a  parabolic  approximation  to  the 
mild  slope  equation  which  accounts  for  combined  refraction  and  diffraction.  The  present 
model  Improves  upon  Yan's  model  in  numerical  efficiency,  by  introducing  structures  and 
most  importantly  by  obtaining  the  radiation  stress  terms  directly  from  the  complex  wave 
amplitude. 
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Figure  1.1:  Computational  domain  of  the  model. 

1.2     General  Description  of  the  Model 

This  section  is  a  description  of  the  numerical  model  that  has  been  constructed  to  calcu- 
late the  wave  and  current  field  within  a  given  domain.  This  domain  usually  is  bordered  by 
a  shoreline  along  one  boundary  and  an  offshore  region  along  the  opposite  boundary.  The 
other  boundaries  are  sufficiently  far  from  the  region  of  interest  (i.e.  a  structure  or  shoal, 
etc.)  so  as  to  not  affect  the  region  of  interest.  This  will  be  discussed  further  in  the  section 
on  boundary  conditions.  The  computational  domain  for  the  purposes  of  this  investigation 
consists  of  a  rectilinear  coordinate  system.  Figure  1.1  shows  the  computational  domain. 

The  model  consists  of  two  parts,  a  circulation  model  and  a  wave  model.  The  circulation 
model  computes  the  mean  water  surface  level  and  the  depth  averaged  mean  currents  using 
the  radiation  stresses  from  the  waves  as  the  forcing  terms  and  includes  bottom  friction, 
advective  acceleration  terms  and  the  lateral  diffusion  of  current  momentum.  The  wave 
model  determines  the  radiation  stress  terms  by  first  solving  for  the  complex  amplitude  as  a 
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function  of  the  total  depth  and  the  depth  averaged  mean  currents,  and  then  computing  the 
radiation  stresses  from  the  complex  amplitudes.  The  model  starts  with  initial  conditions 
of  zero  mean  currents  and  zero  mean  water  surface  elevation.  The  offshore  incident  wave 
height  is  initially  zero  and  then  is  smoothly  increased  up  to  its  given  value  using  a  hyperbolic 
tangent  curve.  This  is  done  so  as  to  minimize  any  transients  due  to  a  shock  start-up. 
The  model  iterates  between  the  wave  model  and  the  circulation  model  until  steady  state 
convergence  is  obtained.  Each  of  the  two  parts  of  the  model  will  be  described  separately. 

The  circulation  model  uses  an  alternating  direction  implicit  (ADI)  method  of  solution. 
Three  equations  are  solved  for  three  unknowns.  The  unknowns  are  fj  the  mean  water  surface 
level,  and  U  and  V,  the  two  components  of  the  depth  averaged  mean  current.  The  equations 
are  the  continuity  equation  and  the  x-  and  y-direction  momentum  equations.  Each  iteration 
of  the  circulation  model  solves  for  the  three  unknowns  on  the  next  time  level.  In  other  words 
each  iteration  advances  the  values  of  fj,  U  and  V  by  a  time  increment,  At.  For  the  purposes 
of  this  investigation  the  model  solves  for  steady  state  conditions  using  the  time  marching 
character  of  the  circulation  model.  The  model  can  be  used  to  solve  for  time  dependent 
wave  fields,  but  this  would  require  the  substitution  of  a  time  dependent  wave  equation  in 
the  wave  model.  The  ADI  method  first  solves,  grid  row  by  grid  row,  for  U  at  the  next  time 
level  and  an  intermediate  value  of  fj,  by  solving  the  x  direction  equation  of  motion  and  the 
continuity  equation  using  an  implicit  numerical  scheme.  This  solution  is  in  terms  of  the 
known  values  of  17, V, and  17  at  the  present  time  level.  Once  all  the  grid  rows  are  solved  in 
the  x  direction,  the  y  equation  of  motion  and  the  continuity  equation  are  solved,  again  by 
an  implicit  scheme,  grid  column  by  grid  column,  for  V  and  fj  at  the  next  time  level  in  terms 
of  U  at  the  new  time  level,  fj  at  the  intermediate  level  and  V  at  the  present  level.  For  each 
of  the  implicit  schemes  in  the  two  directions,  the  governing  equations  yield  a  tridiagonal 
matrix.  The  tridiagonal  matrix  and  the  method  of  solution  will  be  discussed  later. 

The  wave  model  also  yields  a  tridiagonal  matrix  equation.  The  governing  equation  is 
a  parabolic  equation  for  the  complex  amplitude  A.  The  complex  amplitude  contains  phase 
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information  meaning  among  other  things  that  y  variations  of  the  wave  field  are  contained 
within  the  complex  amplitude.  The  term  "parabolic  equation"  indicates  that  there  are 
first  order  derivatives  in  the  x  direction  and  second  order  derivatives  in  the  y  direction. 
The  method  of  solution  is  to  solve  for  the  amplitude  on  the  I-fl  grid  row  in  terms  of  the 
known  amplitude  values  on  the  I  grid  row.  The  solution  scheme  starts  on  the  first  grid  row 
using  the  prescribed  offshore  incident  wave  and  then  advances  grid  row  by  grid  row  to  the 
last  grid  row.  To  advance  to  the  next  grid  row  the  amplitude  values  on  the  1+1  grid  row 
are  implicitly  formulated  using  a  Crank-Nicolson  scheme.  This  will  be  further  discussed  in 
section  4.2  which  describes  the  finite-differencing  scheme  for  the  parabolic  wave  equation. 
Once  the  complex  amplitudes  are  known  the  radiation  stress  terms  are  computed  in  terms 
of  the  absolute  value  of  the  amplitude  and  the  gradients  of  the  amplitude.  This  is  described 
in  section  2.2. 

The  circulation  model  and  the  wave  model  are  the  heart  of  the  program;  however  there 
is  much  more.  In  all  there  are  some  27  subroutines  in  the  program.  These  subroutines  cover 
everything  from  input  and  output  to  convergence  checks  and  flooding  of  the  shoreline.  The 
main  program  is  compartmentalized  into  many  subroutines  so  as  to  facilitate  any  changes 
or  future  upgrading.  For  example,  the  bottom  shear  stress,  a  term  of  importance  in  the 
equations  of  motion,  currently  uses  a  crude  approximation.  If  a  better  approximation  is 
developed  that  is  also  computationally  efficient,  it  will  be  an  easy  task  to  replace  the  present 
subroutine  that  computes  the  bottom  shear  with  the  newer  subroutine.  There  are  many 
places  where  the  model  can  be  improved  with  better  assumptions  and  better  approxima- 
tions. These  will  be  discussed  in  detail  in  the  final  chapter  outlining  recommendations  for 
further  study.  The  modular  construction  of  the  model  also  allows  for  ease  in  substituting 
one  wave  model  for  another. 

Figure  1.2  shows  a  flow  chart  of  the  computer  program.  The  subroutine  INPUT  is  called 
to  initialize  the  values  of  the  unknowns  and  to  input  the  water  depth,  the  offshore  wave 
height  and  direction  and  the  position  of  any  structures.    The  three  subroutines  WAVE, 
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Figure  1.2:  Flow  chart  of  computer  program. 
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SHEAR,  and  LATMIX  compute  the  radiation  stresses,  the  bottom  shear  stress,  and  the 
lateral  mixing  coefficients.  These  terms  are  the  forcing  terms  in  the  equations  of  motion 
and  are  used  in  the  subroutine  CIRC,  which  is  the  circulation  model,  to  solve  for  fj,  U  and  V. 
There  is  a  decision  tree  which  directs  the  program  to  the  three  subroutines  WAVE,  SHEAR, 
and  LATMIX  at  selected  intervals.  This  is  purely  arbitrary  and  is  done  in  the  interest  of 
computer  speed.  It  is  assumed  that  the  forcing  terms  do  not  change  at  a  great  rate  and  can 
thus  be  updated  periodically.  This  updating  is  more  frequent  during  the  start-up  time  while 
the  offshore  incident  wave  height  is  being  increased.  The  CHECK  subroutine  is  to  check  the 
model  for  convergence.  Since  steady  state  conditions  are  assumed,  this  subroutine  indicates 
convergence  when  the  absolute  value  of  the  percentage  difference  between  the  updated  values 
of  U,V,  and  fj  and  the  previous  values  is  less  than  an  arbitrarily  small  specified  value.  The 
UPDATE  subroutine  simply  updates  the  values  of  the  unknowns  for  the  next  iteration.  The 
FLOOD  subroutine  allows  for  the  addition  or  elimination  of  grid  rows  at  the  shoreline  to 
allow  for  beach  flooding  due  to  set-up. 

Figure  1.3  is  a  flow  chart  of  the  circulation  model.  The  subroutines  XCOEF  and 
YCOEF  determine  the  coefficients  of  the  tridiagonal  matrix  equation  for  the  general  cases 
in  the  x  and  y  directions.  The  subroutine  TRIDA  solves  the  tridiagonal  matrix  equation 
using  the  double-elimination  scheme  described  by  Carnahan,  Luther,  and  Wilkes  (1969). 
Subroutines  XCOJ1  and  XCOJN  are  adaptations  of  XCOEF  to  determine  the  coefficients 
of  the  tridiagonal  matrix  for  the  special  case  along  the  lateral  boundary  where  J=l  or 
J=N.  Likewise  YCOIM  is  an  adaptation  of  YCOEF  for  the  I=M  grid  row  which  borders 
the  shoreline.  These  adaptations  are  necessary  because  boundary  conditions  necessitate 
different  numerical  formulation  of  some  of  the  derivative  terms.  Similarly  with  the  presence 
of  an  emerged  impermeable  barrier  bordering  the  shore  side  of  the  JJJ  grid  row,  special 
subroutines  YCJJJ  and  YCJJJ1  are  used  on  the  grid  rows  I=JJJ  and  I=JJJ+1  so  as  to 
properly  include  the  boundary  conditions  for  these  grid  rows,  and  for  shore  perpendicular 
groins  subroutines  XCJGR  and  XCJGR1  incorporate  the  boundary  conditions  for  a  groin 
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Figure  1.3:  Flow  chart  of  circulation  model. 
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Figure  1.4:  Flow  chart  of  wave  model. 

located  between  the  J  =  JGR  and  the  J  =  JGR+1  grid  columns. 

Figure  1.4  is  a  flow  chart  of  the  wave  model.  Subroutine  WAVNUM  determines  the  wave 
number  k,  the  group  velocity  Cg  and  the  intrinsic  frequency  a  at  each  of  the  grid  centers. 
These  terms  are  used  in  the  governing  parabolic  equation.  Subroutine  AMPCO  determines 
the  coefficients  of  the  tridiagonal  matrix  and  CTRIDA  solves  the  matrix  equation  for  the 
complex  amplitude  on  the  next  grid  row.  CTRIDA  is  essentially  the  same  as  TRIDA  except 
that  it  has  been  modified  to  handle  complex  numbers.  The  subroutine  WW  is  called  by  the 
AMPCO  subroutine  in  order  to  determine  the  dissipation  coefficient  following  the  work  of 
Dally  (1980,1987). 

Chapter  2  describes  the  circulation  model  in  greater  detail.  Individual  sections  will 
cover  the  governing  equations  and  the  formulation  of  each  of  the  component  terms.  Chap- 
ter 3  describes  the  wave  model,  the  derivation  of  the  linear  wave  equation,  the  parabolic 
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approximation,  the  manner  of  including  wave  energy  dissipation  due  to  wave  breaking  and 
the  method  of  obtaining  the  radiation  stress  terms  directly  from  the  complex  amplitude. 
The  fourth  chapter  gives  details  of  the  numerical  methods  used  to  obtain  solution  of  the 
governing  equations  and  the  boundary  conditions  employed.  The  fifth  chapter  presents  the 
results  of  applications  of  the  model  and  comparisons  with  other  studies  both  numerical  and 
experimental,  as  well  as  some  analytical  solutions  to  some  idealized  cases.  The  final  chapter 
summarizes  the  study  and  makes  recommendations  for  future  research. 


CHAPTER  2 
CIRCULATION  MODEL 


2.1     Governing  Equations 

The  governing  equations  for  the  circulation  model  are  the  depth  integrated  time  aver- 
aged horizontal  equations  of  momentum 

1    (dS„  ,  dSzv\      ldn 


dV+UdV+VdV  4-  n9f>l       l  l 


and  the  continuity  equation 


1    (dSxy       dSvv\       ldri 


ai  +  r^UD^Ty<yD)  =  0  <23> 


where 


U  =  x  component  of  the  mean  current 

V  =  y  component  of  the  mean  current 

fj  =  mean  water  surface  elevation 

p  =  water  density 

D  =  the  total  water  depth  =  h0  +  fj 

h0  =  local  still-water  depth 

rj  =  the  lateral  shear  stress  due  to  turbulent  mixing 

Tbx  —  x  component  of  the  bottom  shear  stress 
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rby    —    y  component  of  the  bottom  shear  stress 

rtx    =    x  component  of  the  surface  shear  stress 
T»v    =    y  component  of  the  surface  shear  stress 

and  Sxx,  Sxy  and  Svv  are  the  radiation  stress  components  which  arise  from  the  excess 
momentum  flux  due  to  waves.  It  should  be  noted  that  U  and  V  include  the  wave  induced 
mass  transport. 

These  equations  are  obtained  by  integrating  the  local  x  and  y  momentum  equations  and 
the  continuity  equation  over  the  depth  of  the  water  column  and  then  time  averaging  the 
results.  This  is  demonstrated  in  Appendix  A,  which  follows  closely  the  work  of  Ebersole  and 
Dalrymple  (1979).  For  greater  details  of  the  derivation  the  reader  is  referred  to  Ebersole 
and  Dalrymple  (1979)  or  Phillips  (1969). 

These  are  the  same  equations  that  have  been  used  by  other  modelers  such  as  Birke- 
meier  and  Dalrymple  (1976),  Ebersole  and  Dalrymple  (1979),  Vemulakonda  (1984),  and 
Yan  (1987).  Birkemeier  and  Dalrymple  used  an  explicit  numerical  scheme  and  neglected 
the  advective  acceleration  terms  and  the  lateral  mixing.  Ebersole  and  Dalrymple  added  the 
advective  acceleration  and  lateral  mixing  terms  but  still  used  an  explicit  numerical  scheme. 
Vemulakonda  introduced  an  implicit  numerical  formulation  for  these  equations.  Each  of 
these  investigators  used  a  refraction  wave  model  based  upon  the  work  of  Noda  et  al.  (1974) 
which  solves  for  the  wave  height  and  wave  angle  at  given  grid  points.  Yan  (1987)  used  a 
combined  refraction-diffraction  wave  equation  to  determine  the  complex  wave  amplitude 
(magnitude  and  phase)  and  then  solved  for  the  wave  angle  in  terms  of  the  slope  of  the  wa- 
ter surface.  Each  of  these  studies  obtained  the  radiation  stress  components  in  terms  of  the 
wave  height  and  wave  angle.  The  present  study  differs  from  the  previous  ones  in  that  the 
radiation  stress  terms  are  obtained  directly  from  the  complex  amplitude  and  its  gradients. 
This  will  be  discussed  in  section  2.2. 

The  use  of  Eqs.  (2.1-2.3)  involves  many  assumptions.  The  primary  assumption  is 
that  the  flow  field  can  be  represented  in  two  dimensions  using  depth  and  time  averaged 
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equations.  This  assumption  implies  that  vertical  fluctuations  of  the  flow  which  are  zero  in 
the  time  average  can  be  ignored.  Many  details  of  the  flow  may  be  missed  by  this  assumption 
of  purely  horizontal  flow.  Of  course  reducing  the  problem  to  a  two-dimensional  flow  has 
the  distinct  advantage  of  reducing  an  essentially  intractable  problem  to  a  simpler  tractable 
form. 

Some  of  the  terms  in  Eqs.  (2.1-2.3)  also  involve  certain  assumptions.  These  terms 
are  the  bottom  friction,  the  surface  wind  stress  and  the  lateral  mixing.  These  are  all 
phenomenon  that  are  not  completely  understood  so  that  assumptions  are  made  concerning 
the  complex  forces  that  are  represented.  These  assumptions  usually  involve  the  use  of 
empirical  coefficients.  In  this  study  the  wind  forces  are  ignored  so  the  assumptions  involved 
with  the  surface  stress  will  not  be  discussed  in  detail  here.  It  would  be  a  relatively  simple 
task,  though,  to  add  a  wind  surface  stress  to  the  model,  using  a  quadratic  shear  stress  as 
formulated  by  Birkemeier  and  Dalrymple  (1976)  and  also  used  by  Ebersole  and  Dalrymple 
(1979)  following  the  work  of  Van  Dorn  (1953)  in  determining  the  value  of  the  empirical 
coefficient. 

The  interaction  of  an  oscillatory  fluid  flow  with  a  real  bottom  is  very  complex.  Many 
physical  phenomena  may  be  at  work.  With  sand  bottoms  there  can  be  percolation  into  the 
sand  and  frictional  effects  resulting  in  wave  energy  dissipation.  With  mud  bottoms  where 
the  mud  is  essentially  a  dense  fluid  layer  there  can  be  damping  and  internal  waves  along 
the  mud-water  interface.  Wave  energy  losses  due  to  bottom  effects  are  not  included  in  the 
model  but  can  easily  be  incorporated  into  the  wave  dissipation  subroutine.  Resuspension 
and  settling  of  particles  also  constitutes  a  complex  phenomenon  with  density  changes  and 
momentum  transfers.  All  of  these  effects  are  glossed  over  in  the  model  by  assuming  a  rigid 
impermeable  bottom.  As  with  all  of  the  models  discussed  above  the  bottom  friction  term  is 
represented  in  the  quadratic  form.  Noda  et  al.  (1974),  Birkemeier  and  Dalrymple  (1976)  and 
Vemulakonda  (1984)  all  used  the  weak  current  formulation  as  proposed  by  Longuet-Higgins 
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(1970a)  in  which 

Tbx    =    pF  \uorb\  U  (2.4) 

Thy     =    pF  \uorb\  V  (2.5) 

where  F  is  a  drag  coefficient  (of  the  order  of  .01)  and  |uorj,|  is  the  time  average  over  a  wave 
period  of  the  absolute  value  of  the  wave  orbital  velocity  at  the  bottom.  From  linear  wave 
theory  it  is  found  that 

where  H  is  the  wave  height  and  T  is  the  wave  period,  k  is  the  wavenumber,  and  h  =  D  is 
the  total  water  depth.  Ebersole  and  Dalrymple  used  a  more  complete  quadratic  formulation 


Tb  =  pFut\ut\  (2.7) 

where  the  overbar  represents  the  time  average  over  a  wave  period  and  ut  is  comprised  of  both 
the  mean  current  and  the  wave  orbital  velocities.  The  total  velocity  vector  uj  is  expressed 
as 

tTt  =  (U  +  u  cos  9)t  +  (V  +  u  sin  0)j  (2.8) 

so  that  its  magnitude  is  given  by 

|tT,|  =  VU2  +  V2  +  ti2  +  2Uu  cos  9  +  2Vu  sin  0  (2.9) 

The  wave  orbital  velocity  u  is  expressed  as 

u  =  tim  cos  at  (2.10) 

where  tim  is  the  maximum  wave  orbital  velocity  at  the  bottom  which  is  found  to  be 

(2.11) 


Um  = 


T  sinh  A:Ji 
The  x  and  y  components  of  the  time  averaged  bottom  friction  are  thus 


pF   f2x 
Hz    =     ^y     {U  +  um  cosfl  cos  at)  -\ut\d{at)  (2.12) 

pF  f2* 
nv     =     ^J      (V  +  umsin0cos at)  ■  \ut\d{at)  (2.13) 
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where  \ut\  is  written  as 


l«3|  =  \/^2  +  V2  +  u^  cos2  at  +  2£7um  cos  at  cos  0  +  2Kum  cos  at  sin  0 


(2.14) 


The  integrals  in  Eqs.  (2.12)  and  (2.13)  are  then  evaluated  by  a  16  point  summation  using 
Simpson's  Rule. 

In  the  present  study  several  different  formulations  are  used  and  compared.  The  first  is 
the  formulation  given  by  Longuet-  Higgins  (1970a)  and  used  by  several  previous  investiga- 
tors. The  second  is  a  modification  of  the  Longuet-Higgins  formulation  with  partial  inclusion 
of  current  strength  in  the  friction  formula. 

rbx    =    pF  |uori|  U  +  pF\U\U  (2.15) 

rbv    =    pF  \uorb\  V  +  pF\U\V  (2.16) 

It  is  thought  that  this  formulation  has  some  value  since  in  the  lee  of  structures  where  there 
is  relatively  weak  wave  activity,  there  may  still  be  relatively  strong  currents  generated  by 
the  gradients  of  the  set-up  outside  of  the  sheltered  region.  The  third  formulation  is  similar 
to  that  of  Ebersole  and  Dalrymple  in  that  the  full  quadratic  expression  is  used.  However  the 
wave  orbital  velocities  at  the  bottom  are  expressed  as  gradients  of  the  complex  amplitude 
and  the  resulting  integral  is  determined  using  an  eight  point  Gaussian  integration.  It  i 
found  that  an  eight  point  Gaussian  integration  of  a  cosine  squared  integration  over  an 
interval  of  2tt  is  accurate  to  within  1  percent. 

The  quadratic  formulation  of  the  bottom  friction  is  used 


is 


Tb  =  pFut\ut\  (2.17) 

where  ut  is  composed  of  both  the  mean  current  and  the  wave  orbital  velocities.  The  wave 
orbital  velocities  are  expressed  as  the  real  part  of  the  gradient  of  the  complex  velocity 
potential  $  which  is  expressed  as 

g=      A(X,y)coshk(h+z)e_iwt 

a  cosh  kh  K  '     > 
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where  A  is  the  complex  amplitude,  u  is  the  absolute  frequency,  a  is  the  intrinsic  or  relative 
frequency  following  the  current,  and  h0,  h  and  *  are  as  previously  explained.  Equations 
governing  $  and  a  parabolic  approximation  solution  for  A  are  presented  in  the  next  chapter. 
The  total  velocity  vector  ut  is  expressed  as 

•-[■♦■(SM'+Mie)]*  m 

where  5R  indicates  the  real  part,  so  that  its  magnitude  is  given  by 


«'-r+v,+ra(S)+»*(S)+KS)],+K5)r  <-» 

The  x  and  y  components  of  the  time  averaged  bottom  friction  are  thus 

*    =    pFfJo[V  +  *{i%)]-Wdt  (222) 

As  in  Ebersole  and  Dahymple,  the  assumption  is  made  that  only  the  gradients  of  A  are 
retained  and  the  gradients  of  the  depth  and  spatially  varying  quantities  such  as  k  and  a 
are  ignored.  This  is  essentially  the  assumption  of  a  flat  bottom  and  very  slowly  varying 
currents. 

Figure  2.1  plots  the  magnitude  of  the  longshore  current  as  a  function  of  the  distance 
offshore  for  the  various  friction  formulations  discussed  above.  These  velocity  profiles  were 
obtained  numerically  for  an  eleven  second  wave  with  a  deep  water  wave  height  of  1.028 
meters  (1  meter  at  offshore  grid  of  the  computational  domain,  using  a  AX  of  10  meters 
and  30  grid  rows)  and  a  deep  water  wave  angle  of  17.28  degrees  (10  degrees  at  the  model 
boundary)  approaching  a  1:25  plane  beach.  A  friction  factor,  F,  of  .01  was  used  for  each  case. 
Using  a  friction  formulation  that  ignores  the  contribution  of  the  waves  results  in  a  very  weak 
friction  and  thus  a  large  longshore  current.  Using  the  Longuet-  Higgins  formulation  gives  a 
friction  slightly  less  than  the  more  proper  time  averaged  formulation,  while  a  combination  of 
the  Longuet- Higgins  formulation  and  the  strong  current  model  (no  wave  contribution)  yields 
a  slightly  stronger  friction.  For  practical  numerical  considerations  the  latter  formulation  is 
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Figure  2.1:  Longshore  current  profiles  for  different  friction  formulations  for  an  11  second, 
1.028  meter  wave  at  17.28  degrees  to  a  1:25  plane  beach. 

most  often  used  since  it  requires  less  computation  time  and  more  importantly  the  Gaussian 
integration  involving  the  real  part  of  the  gradient  of  complex  amplitudes  at  times  produces 
stability  problems. 

The  lateral  mixing  term  in  Eqs.  (2.1)  and  (2.2)  also  involve  several  simplifying  assump- 
tions. The  first  assumption  is  found  in  the  derivation  of  the  two  depth  integrated  equations 
of  motion  where  it  is  assumed  that  the  turbulent  shear  stress  term  -puV  is  independent 
of  depth.  The  other  assumption  is  that  the  process  of  momentum  transfer  due  to  turbulent 
fluctuations  can  be  represented  by  a  product  of  a  mixing  length  parameter  and  derivatives 
of  the  mean  current.  This  is  what  is  done  in  the  present  study  and  also  in  previous  studies. 
Ebersole  and  Dalrymple  (1979),  Vemulakonda  (1984),  and  Yan  (1987)  all  used  the  formu- 
lation given  by  Longuet-Higgins  (1970b)  in  which  he  developed  an  analytical  expression  for 
the  cross  shore  distribution  of  the  longshore  current  velocities  using  a  linearized  friction 
formulation  as  described  above,  linear  waves  at  an  angle  to  the  normal  of  a  plane  beach, 
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and  a  spilling  breaker  assumption.  Longuet-Higgins  wrote 

dU  BV\ 


in  which  ex  and  e„  are  mixing  length  coefficients. 

Longuet-Higgins  reasoned  that  the  mixing  length  parameter  ex  must  tend  to  zero  as 
the  shoreline  was  approached.  He  assumed  that  ex  is  proportional  to  the  distance  from  the 
shoreline,  X,  multiplied  by  a  velocity  scale  which  he  chose  as  the  shallow  water  celerity 
\fgk.  Thus  ex  is  written  as 

ex  =  NXy/gh  (2.24) 

where  N  is  a  dimensionless  constant  for  which  Longuet-Higgins  chose  the  limits  as 

0   <    N   <   0.016  (2.25) 

Obviously,  for  physical  reasons  there  needs  to  be  some  limit  on  the  scale  of  these  eddies. 
Ebersole  and  Dalrymple  arbitrarily  chose  the  value  of  ex  near  the  breaker  line  to  be  the 
maximum  value  and  for  ex  to  be  constant  offshore  from  the  breaker  line.  The  same  limit 
upon  ex  is  used  in  this  study.  The  coefficient  e„  is  assumed  to  be  a  constant  and  it  is 
arbitrarily  assigned  the  maximum  value  of  ex. 

In  the  present  study  the  further  assumption  is  made  in  computing  rj  that  the  cross 
derivatives  are  negligible  (a  common  assumption  in  fluid  mechanics).  Therefore,  the  lateral 

mixing  terms  become 

1  dn  d   (    au\ 

^  =  "^  PW  (226) 

in  the  x  direction  and 

I  dn  d    (     dV\ 

-pJ-x^-di^)  (227) 

in  the  y  direction.  This  means  that  the  lateral  mixing  term  is  represented  by  a  purely  lateral 
diffusion  of  momentum.  That  the  lateral  mixing  term  is  a  lateral  diffusion  term  has  some 
significance  in  that  the  diffusion  term  is  expressed  explicitly  and  thus  there  results  a  time 
step  limitation  to  an  implicit  model.  This  will  be  discussed  in  section  4.1  dealing  with  the 
finite  differencing  of  the  equations  governing  the  circulation  model. 
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2.2     Radiation  Stresses 

In  a  series  of  papers  in  the  early  1960s  Longuet-Higgins  and  Stewart  (1962,  1963,  1964) 
introduced  and  developed  the  mathematics  of  the  radiation  stresses.  These  stresses  which 
can  be  conceptually  thought  of  as  the  "flux  of  excess  momentum"  can  be  obtained  through 
integration  of  the  momentum  equations  over  the  water  column  as  is  shown  in  Appendix  A. 
Using  linear  wave  theory  the  radiation  stress  components  are  found  to  be 


Szz    —    E 

Oyy         =         E 


n(cos20+l)-l]  (2.28) 

n(sin20  +  l)-l]  (2.29) 


Sxy    =    Encos6sai9  (2.30) 

where  E  is  the  wave  energy  equal  to  \ pgH2,  n  is  the  ratio  of  the  group  velocity  to  the  wave 
celerity  and  0  is  the  angle  the  wave  rays  make  with  the  x  axis. 

If  the  wave  height  and  wave  angle  are  known  at  each  grid  point,  the  radiation  stresses 
and  their  gradients  are  readily  determined  for  use  in  the  governing  equations.  To  determine 
the  wave  height  and  angle  at  each  grid  point  Birkemeier  and  Dalrymple  (1976),  Ebersole 
and  Dalrymple  (1979)  and  Vemulakonda  (1984)  each  used  the  following  scheme  which  is 
based  upon  the  irrotationality  of  the  wave  number  and  an  equation  for  the  conservation  of 
wave  energy.  The  irrotationality  of  the  wave  number  determines  changes  in  the  wave  angle 
due  to  refraction  while  the  energy  equation  determines  changes  in  the  wave  height  due  to 
shoaling  and  refraction.  This  scheme  was  first  introduced  by  Noda  et  al.  (1974). 

The  irrotationality  of  the  wave  number 

V   x   £=0  (2.31) 

where 

*  =  |fc|  cos  0i+  \k\  sin  0]  (2.32) 

is  expanded  in  Cartesian  coordinates  to  obtain 

ode   ,        »d6       ■    «  1  d\k\  1  d\k\ 


25 
For  waves  on  a  current  the  dispersion  relation  is 


u 


=  a  +  k-U  =  y/g\k\  tanh  \k\h  +  \k\  cos  0  +  \k\  sin  0  (2.34) 

Equation  (2.34)  is  differentiated  to  eliminate  ^^1  and  jfj^gl  from  Eq.  (2.33).  Using  a 
forward  difference  derivative  in  x  and  a  backward  difference  derivative  in  y,  $tJ  is  obtained 
as  a  function  of  cos  0 ,  sin  0,  U,  V,  k,h  and  their  derivatives  as  well  as  0  at  the  adjoining  grids. 
The  sin  0  and  cos  0  terms  are  expressed  in  terms  of  the  surrounding  four  grids  using  a  2nd 
order  Taylor  expansion.  The  solution  is  an  iterative  process  alternating  solution  of  the  0 
equation  with  solution  of  the  dispersion  relationship  (2.34)  for  k  =  \k\  which  is  accomplished 
by  a  Newton-Raphson  method  until  acceptable  convergence  is  obtained  for  k  and  0. 

Next  an  energy  equation  is  solved  for  the  wave  height.  A  conservation  of  energy  equation 
is  given  by  Phillips  (1969  ed.  Eq.  3.6.21)  as 

£+£lHR  +  <U}  +  *«i.O  (2.35) 

Assuming  steady  state  conditions  and  expanding  in  Cartesian  coordinates  gives 

where 

*«.     =      -gS„        =     n(cos2  0  +  1)-  - 
svv     =     -^svv        =    "(sin2  0  +  1)  -  - 


fXy  —  tVz  —  -p  Sxy  =  —  Sxy    =    ncos0sin0 


Defining 


where 


QiJ  =  \{'k{U  +  C'coa^  +  ^V  +  C'  8in  9)  +  *}  (2-37) 


a  = 


dU_  dU  dV  _     dV] 


(2.38) 
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and  noting  that  U,V,Cg,  and  9  are  known  quantities;  taking  forward  difference  derivatives 
in  x  and  backward  difference  derivatives  in  y  and  solving  for  the  wave  height,  Htj,  at  the 
i,j  grid  yields 

H'J  =  — v+c.  rin< — v+c.  co,  e*X„  (2-39) 

The  computation  is  a  row  by  row  iteration  proceeding  from  large  i  (offshore  in  the 
coordinate  system  of  the  above  cited  authors)  to  small  i.  During  the  computation  of  each 
Hi j  the  breaking  height  is  also  calculated  and  if  H, \j  exceeds  this  value,  then  the  breaking 
wave  height  is  used.  Several  breaking  height  indices  are  possible.  Miche  (1944)  gave  as  the 
theoretical  limiting  wave  steepness  the  condition 

XI 

-^  =  .142tanh(fct/i6)  (2.40) 

where  the  b  subscript  indicates  breaking  conditions.  Le  Mehaute  and  Koh  (1967)  indicated 
from  experimental  data  that  the  limiting  steepness  is 

-j4  =  .12tanh(M*)  (2.41) 

or 

Hb  =  ^-.12t&nh[kbhb)  (2.42) 

This  is  the  breaking  height  criterion  used  by  Noda  et  al.    (1974).    Perhaps  the  simplest 
breaking  height  criterion  is  the  spilling  breaker  assumption 

#4  =  -78/it  (2.43) 

which  is  used  by  Vemulakonda  (1984)  and  Vemulakonda  et  al.  (1985). 

The  scheme  described  above  of  determining  the  wave  angle  6  from  the  irrotationality 
of  the  wave  number  (2.33)  and  the  dispersion  relation  (2.34)  and  the  wave  height  from 
the  energy  equation  (2.35)  has  its  limitations.  Principally  this  scheme  does  not  take  into 
account  the  effects  of  diffraction.  This  limitation  is  surmounted  by  use  of  a  parabolic  wave 
equation.  However  this  presents  another  problem  in  that  the  wave  angle  is  not  determined 
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in  a  parabolic  equation  solution.  Only  the  complex  amplitude,  which  contains  the  mag- 
nitude and  the  phase  of  the  wave,  is  obtained.  Yan's  (1987)  solution  to  this  problem  is 
to  obtain  0  directly  from  the  complex  amplitude  by  assuming  that  the  slope  of  the  water 
surface  resulting  from  the  waves  indicates  the  direction  of  wave  propagation.  This  value 
for  0  and  twice  the  absolute  value  of  the  complex  amplitude  for  the  wave  height  are  then 
used  in  Eqs.  (2.28-2.30)  to  obtain  the  radiation  stresses.  The  value  thus  obtained  for  6  is 
then  also  used  in  the  dispersion  relationship.  For  plane  waves  this  works  fine;  however,  in 
diffraction  zones  where  short-crested  wave  are  encountered  problems  may  arise,  since  the 
instantaneous  slope  of  the  wave  at  a  particular  location  can  be  other  than  in  the  direction 
of  wave  propagation. 

In  this  report  the  radiation  stresses  are  obtained  directly  from  the  complex  amplitude 
and  its  gradients  using  equations  developed  by  Mei  (1972)  and  repeated  in  his  book  (1982). 
The  basis  of  the  derivation  is  to  use  the  definition  of  the  velocity  potential  to  substitute  for 
the  wave  induced  velocities  in  the  definition  for  the  radiation  stresses. 

For  an  inviscid  incompressible  fluid  the  instantaneous  equations  of  motion  are 

P(d*+''      7     =     -VP-P9e*  (2.44) 

V-f    =    0  (2.45) 

where  q  represents  the  velocity  components,  p  is  the  pressure  and  ex  is  the  z  or  vertical  unit 
vector.  Denoting  the  horizontal  components  of  fas  ua,  a  =  1, 2  and  the  vertical  component 
of  q  as  w,  the  boundary  conditions  are 

P    =    0  (2.46) 

dn  dn 

on  z  =  n,  where  n  is  the  instantaneous  free  surface,  and  the  bottom  boundary  condition 

dh 

q  is  assumed  to  be  composed  of  a  time  averaged  part  q  and  a  fluctuating  part  q'\  i.e. 

?=${x,V,z)  +  q'{x,V>z,t)  (2.49) 
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Integrating  a  horizontal  component  of  (2.44)  vertically  over  the  depth  of  the  water  column 
from  z  =  -h  to  z  =  17,  using  the  boundary  conditions  (2.46-2.48)  and  then  taking  time 
averages  (the  same  procedure  as  in  appendix  A)  the  following  definition  of  the  radiation 
stresses  is  obtained: 


Sap  =  pj°_h  u'au'0dz  +  6af)  \PhPdz  -  f(fj  +  h)>  +  ^{^ 


(2.50) 


where  the  overbar  is  used  to  indicate  the  time  averaged  values  and  6afi  is  the  Kronicker 
delta  function 


f  1  for  a  =  3 

"   =  I  ° 


'f  ~  \  0  for  a  *  8  (251) 

Mei  (1972)  was  able  to  express  Eq.  (2.50)  in  terms  of  the  general  complex  amplitude 
function  B  by  relating  it  to  the  velocity  potential  and  the  free  surface  height  by 

#  -  -ft*=ftife+d*-  (2.52) 


r,    =     Be~iut 


a  cosh  kh 

(2.53) 


An  expression  for  p  is  obtained  by  integrating  the  vertical  component  of  Eq.  (2.44) 

dw  dw  dw  dp 

PH+pU(,d^  +  pW-dz-  =  -Tz-pg  <2-54) 

Integrating  from  any  depth  z  to  the  surface,  using  Leibnitz's  rule  for  interchanging  the 
order  of  differentiation  and  integration,  employing  the  two  surface  boundary  conditions, 
Eqs.  (2.46)  and  (2.47)  and  the  continuity  Eq.  (2.45),  Eq.  (2.54)  yields: 

p(x,  y,  z,  t)  =  pg(n  -z)  +  p—j\dz  +  p^-  f"  Upw  dz  -  pw\  (2.55) 

Taking  time  averages  gives 


P{x,y,z)  =  pg(f)-z)  +  pJ^    —JL— dz  -  p(w'M)*  (2.56) 

Substitution  of  Eqs.  (2.52),(2.53)  and  (2.56)  into  (2.50)  and  making  the  assumption  that 
horizontal  derivatives  of  k,h  and  a  can  be  neglected  (this  is  essentially  a  mild  slope  as- 
sumption), the  following  is  obtained  after  straightforward  but  tedious  calculations  which 
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are  shown  in  Appendix  B: 

*»  [W^kh  +  ^^  -  ^B\2)^(2khcotb2kh  -  1)]}  (2.57) 

where  B*  denotes  the  complex  conjugate  of  B.  In  obtaining  Eq.  (2.57)  use  was  made  of  the 
flat  bottom  Helmholtz  equation 

V2hB  +  k2B  =  0  (2.58) 

The  general  complex  amplitude  function  B  is  related  to  the  complex  amplitude  A  which 
is  solved  for  in  the  wave  model  by 


B  =  Ae^flco,Sdx) 


(2.59) 


There  are  some  problems  with  Eq.  (2.57).  One  may  question  the  assumption  of  a 
flat  bottom;  however  the  inclusion  of  the  slope  and  gradients  of  the  wave  number  and  the 
intrinsic  frequency  leads  to  intractable  mathematics.  The  use  of  the  localized  flat  bottom 
is  also  implicit  in  the  derivation  of  Longuet-Higgins  and  Stewart  to  obtain  Eqs.  (2.28-  2.30) 
since  in  carrying  out  the  integrations  involved  it  is  assumed  that  horizontal  derivatives  of 
the  velocity  potential  are  limited  to  derivatives  of  the  phase  function.  For  plane  waves 
using  a  spilling  breaker  assumption,  results  using  Eq.  (2.57)  differ  little  from  those  using 
Eqs.  (2.28-2.30);  however  in  diffraction  zones  Eq.  (2.57)  is  superior  for  reasons  mentioned 
previously. 

The  major  problem  with  the  use  of  Eq.  (2.57)  is  also  a  major  problem  associated  with 
the  use  of  Longuet-Higgins's  formulation.  This  problem  is  that  in  both  formulations  the 
initiation  of  wave  breaking  implies  the  dissipation  of  energy  and  thus  the  forcing  of  set- 
up and  longshore  currents  starting  at  the  breaker  line.  The  physical  reality  is  somewhat 
different  in  that  in  the  initial  stages  of  breaking,  energy  is  not  immediately  dissipated  but 
rather  is  transformed  from  organized  wave  energy  to  turbulent  energy  and  to  a  surface  roller 
(Svendsen  1984).  There  is  a  space  lag  between  the  initiation  of  breaking  and  the  dissipation 
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Figure  2.2:  Experimental  results  for  wave  set-up  and  set-down.  Reprinted  from  Bowen,  In- 
man,  and  Simmons  (1968),  Journal  of  Geophysical  Research,  vol.  73(8)  page  2573,  copyright 
by  the  American  Geophysical  Union. 

of  energy.  This  can  be  seen  in  figure  2.2  which  gives  the  results  for  the  experiments  of  Bowen 
et  al.  (1968).  Analytical  and  numerical  solutions  for  the  mean  water  surface  elevation  will 
give  a  set-up  starting  at  the  break  point  while  the  solid  line  in  the  lower  graph  shows 
that  set-up  measured  in  the  experiment  begins  shoreward  of  the  break  point.  This  will  be 
discussed  further  in  Chapter  5. 


CHAPTER  3 
WAVE  MODEL 


3.1     Governing  Wave  Equation 


Several  methods  have  been  used  to  derive  the  governing  equation  for  the  propagation 
of  waves  over  varying  topography  and  varying  currents.  Among  the  methods  used  are  per- 
turbation techniques,  multiple  scales,  a  Lagrangian  and  the  use  of  Green's  second  identity. 
With  the  assumption  of  a  mild  slope  and  weak  current  conditions  each  of  these  methods 
should  yield  the  same  equation.  For  illustrative  purposes  the  method  of  using  Green's  sec- 
ond identity  will  be  used  here  to  derive  a  linear  equation.  A  Lagrangian  derivation  of  the 
same  equation  is  presented  in  Appendix  C.  The  linear  equation  is  used  since  second  order 
waves  exhibit  singularities  as  the  water  depth  goes  to  zero. 

A  velocity  potential  <j>  is  assumed  such  that  the  water  particle  velocities  are  given  by 
V^,  where  V  is  the  three-dimensional  gradient  operator 

_        d  - .       d  -.       d  t 

v  =  dit  +  TyJ  +  Tzk  <M) 

and  i,j,  and  k  are  the  unit  vectors  in  the  x,y,  and  z  directions,  respectively.  The  velocities 
are  a  superposition  of  steady  currents  and  wave  induced  motion.  The  velocity  potential  is 
given  by 

4>  =  4>o  +  tjl  (32) 

where  <t>0  is  defined  such  that  its  gradient  yields  U,V,  and  W  which  are  the  steady  current 
terms,  and  tf  is  the  wave  component  of  the  velocity  potential.  An  undefined  scaling  e 
is  used  so  as  to  separate  the  wave  and  steady  current  part  of  the  velocity  potential.  The 
wave  part  of  the  velocity  potential  is  composed  of  two  parts:  /,  which  contains  the  depth 
dependency  and  £,  which  is  a  function  of  the  horizontal  coordinates  and  time.  The  depth 
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dependency  /  is  given  by  linear  theory  as 

.      cosh k(h0  +  z) 

1  "        cosh  kh  (S'S) 

where  h0  is  the  depth  of  the  bottom  referenced  to  the  still  water  line,  and  h  is  the  total 

water  depth  and  k  is  the  wave  number  vector.  It  should  be  noted  that  /  is  also  a  function 

of  the  horizontal  coordinates  since  A;  and  h0  vary  with  horizontal  position. 

For  incompressible,  irrotational  flow  the  governing  equation  for  the  velocity  potential 

is  the  Laplace  equation 

VV  =  0  h0  <  z  <  v  (3.4) 

where  r/  is  the  instantaneous  water  surface  elevation. 
Using  a  horizontal  gradient  operator  V*  defined  by 

the  Laplace  equation  is  written  as 

V\<t>  +  <t>xz  =  0  h0<z<ri  (3.6) 

The  governing  equation  is  complemented  by  the  following  boundary  conditions.  On  the  free 

surface  the  kinematic  free  surface  boundary  condition  is  given  by 

dr] 

jj  +  Vh<t>  •  Vfctj  =  <j>,  z  =  n  (37) 

The  dynamic  free  surface  boundary  condition  with  zero  atmospheric  pressure 

&+-(V^)2  +  ff*  =  0  z  =  r,  (3.8) 

and  the  bottom  boundary  condition 

-  ^h<f> '  V/,/»0  =  <f>x  z=-h0  (3.9) 

The  derivation  technique  is  to  write  Green's  second  identity  for  ^  and  /  as  defined  by 
Eqs.  (3.2)  and  (3.3)  respectively. 

J_     f<t>**dz  -  j     +f„  dz  =  [ft,  -  +ft]lho  (3.10) 
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Then  Eq.  (3.6)  is  used  to  substitute  for  <f>Mt  in  the  first  term  of  Eq.  (3.10)  and  the  boundary 
conditions  are  used  to  substitute  for  <f>„  at  the  mean  free  surface  and  at  the  bottom  in  the 
third  term  of  Green's  identity.  Equation  (3.2)  is  used  and  terms  of  order  e  only  are  retained 
in  order  to  obtain  the  governing  equation  for  4>. 
In  order  to  obtain  <$>t  at  fj  it  is  assumed  that 


r]  =  T]  +  efj. 


(3.11) 


The  kinematic  free  surface  boundary  condition,  Eq.  (3.7),  is  then  expanded  in  a  Taylor 
series  about  z  =  fj, 


^j  +  V^-V^ 


„3 


dz 


Ji  +  Vh<t>-  V/,f?  -  <f>x 


=  0 


M=ff 


(3.12) 


solving  for  <f>z  |^  using  Eqs.  3.2  and  3.11  while  retaining  only  terms  of  order  e  yields 


A  total  horizontal  derivative 


-£  +  U-VhfJ  +  Vh4>-  Vkij  -  rjW» 


D        d 


*=l 


(3.13) 


(3.14) 


is  introduced  so  that  the  first  two  terms  on  the  right  hand  side  of  Eq.  (3.13)  can  be  written 
as  the  total  horizontal  derivative  of  fj.  The  final  manipulation  of  Eq.  (3.13)  is  to  recognize 
that  -^  at  the  surface  is  the  horizontal  divergence  of  the  horizontal  current  (=  Vh'U). 
Thus 


Dfj 
<t>*  \fi=  -ft  +  Vh<f>  •  Vhfj  +  rjVk-U 


(3.15) 


fj  is  eliminated  from  Eq.  (3.15)  by  expanding  the  dynamic  free  surface  boundary  condition, 
Eq.  (3.8),  in  a  Taylor  series  about  z  =  fj. 


h  +  ^tf  +  gz 


=    91  + 


U  +  |(V*)2' 

Using  Eq.  3.2  to  substitute  for  <f>  the  order  t  terms  yield  a  solution  for  fj 

9  L  J  9  Dt 


(3.16) 


(3.17) 
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Substituting  into  Eq.  (3.15)  gives 

+'  U=  —gmi  ~  jF»-U)^  +  V»+-Vk<l  (3.18) 

The  first  term  of  the  Green's  identity  using  Eq.  (3.6)  becomes 

j_ho  f*»dz  =  -  P_K  fVltdz  =  -  jT    fVh  •  (Vk<f>)dz  (3.19) 

Now 

Vh4>=0  +  efVhj>  +  e$Vhf  (3.20) 

and 

Va  •  (Vk<f>)  =  Vh  •  ff  +  «V*  •  (/  VA£)  +  eVh£  •  V*/  +  efalf  (3.21) 

Retaining  only  terms  of  order  t  and  dropping  the  last  term  since  it  is  second  order  in  the 
slope  the  integral  is  equal  to 

J-hB  fVh'V  Vh^dz  ~  j_h  f  V"*  V*/  dz  (3-22) 

These  two  terms  are  handled  using  Leibnitz's  rule  of  integration. 

Va/      f{fVk4)dz    =     f    Vkf{fVk$)dz+  ["    fVh{fVh4>)dz 

°  J -ho  J -h0 

+V^/2  |„  Vfc£  +  Vhh0f*  Uo  Vk$  (3.23) 


Using  this  relationship  and  recognizing  that  /  |„  is  equal  to  one  and  the  integral  over  the 
depth  of  f2  is  CCt/g,  where  C  is  the  wave  celerity,  Cg  is  the  wave  group  velocity  and  g  is 
the  gravitational  constant,  the  first  term  is  equal  to 

CC 
~  V*  '  (-7* V^)  +  Vhfj  Vj  +  Vhh0f*  \.ko  Vh$  (3.24) 

where  the  e  has  been  dropped. 

The  second  term  of  Eq.  (3.10)  is  evaluated  quite  easily. 

m 


-  f     <t>  fzzdz  =  -  f     4>0fttdz  -  €  F    ft  fzzdz 

J~h°  J-ho  J-h0 


(3.25) 
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The  first  integral  above  is  dropped  since  there  are  only  steady  current  terms.  The  remaining 
term,  dropping  the  e,  is 

=  -  f    Si S.,dz  =  -k2$  I"    S2dz  =  -k2™±j>  (3.26) 

The  third  term  in  Eq.  (3.10)  is 

St,  \lh  =  S4>.  U  S**  \-K  (3.27) 

The  first  term  on  the  right  hand  side  is  equal  to  <f>t  |„  since  /  [,=  1  and  has  already  been 
determined.  The  second  term  on  the  right  hand  side  of  Eq.  (3.27)  is  easily  determined  using 
the  bottom  boundary  condition  Eq.  (3.9) 

-  f+m  \-h0=  fVh<t>  •  Vhh0  \.ko  (3.28) 

Using  Eq.  (3.20)  this  becomes 

-  S4m  \-h0=  fU  •  Vhh0  \-K  +ef2  \.ho  Vhj>  ■  Vhh0  +  ef$Vhf  .  Vhh0  \.ho  (3.29) 

The  first  term  is  dropped  since  it  not  order  e  and  the  third  term  is  clearly  higher  order  in 
the  slope.  Thus,  again  dropping  the  e 

-  S<t>z  \-h0=  f2  \-h0  Vh$  ■  Vhh0  (3.30) 

This  term  cancels  the  bottom  boundary  term  resulting  from  the  Leibnitz  integration  above. 
The  fourth  term  of  Eq.  (3.10)  is  also  easily  determined. 

a  t   ifl  i  i  sinh  k(h0  +  z)  ,-  a2 

~  +f.  £*.-  -#      Jh°kh     l  ]<_,=  -#tanh  M  |„=  -°-4>  |,  (3.31) 

Evaluating  <f>  at  the  mean  water  surface  and  retaining  only  terms  of  order  e  yields 

a2  - 
-<t>f*\\0=-—<i>  (3.32) 

Collecting  terms  from  (3.24) ,(3.26) ,(3.18) ,(3.30)  and(3.32)  yields  the  following  linear 
governing  equation  for  <j>. 

D26  m  D6 

DF  +  {Vh'  ^JTt  ~  V*(CCi**+)  +  (**  ~  #CC,)i  =  0  (3.33) 

This  equation  was  first  derived  by  Kirby  (1984). 
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3.2     A  Parabolic  Approximation 

Equation  (3.33)  is  a  hyperbolic  wave  equation.  It  can  be  altered  to  an  elliptic  equation 
by  assuming  that  any  time  dependency  is  solely  in  the  wave  phase.  Before  doing  this  an 
energy  dissipation  term  DD  is  introduced  as  follows.  It  is  assumed  that  Eq.  (3.33)  holds 
for  the  case  of  no  energy  dissipation;  then  for  the  case  with  dissipation  the  left  hand  side  is 
balanced  by  any  energy  that  is  dissipated.  The  governing  equation  is  then  written  as 

D?  +  (V"  '  V)Wt  "  V*(CC«V^)  +  (**  -  #CCg)4>  =  DD  (3.34) 

In  the  derivation  of  the  governing  hyperbolic  equation  using  a  Lagrangian  as  presented  in 
Appendix  C  the  DD  term  can  be  viewed  as  a  virtual  work  term. 

Following  Kirby  (1983)  the  DD  term  is  assumed  to  be  of  the  form 

D6 

DD  =  ~W^7  (3-35) 

where  w  is  an  undefined  coefficient  indicating  the  strength  of  the  dissipation.  Eventually 
the  w  term  will  be  related  to  the  energy  dissipation  due  to  wave  breaking  following  the  work 
of  Dally,  Dean  and  Dalrymple  (1984). 

Using  Eq.  (3.35)  in  Eq.  (3.34)  and  making  the  assumption  that  the  only  time  dependency 

is  in  the  phase,  i.e. 

dd> 

-Qi  =  -""*  (3.86) 

reduces  the  hyperbolic  equation  (3.34)  to  the  elliptic  form 

-2tW  •  Vrf  +  U  •  Vh(U  ■  Vh$)  +  (V„  •  &)(&  •  Vj)  -  V„  •  (CCgVj) 

+{<r2  -  w2  -  k2CC„  -  iu>(Vfc  •  U)}j>  =  iaw4>  (3.37) 

where  only  the  phase  contribution  to  the  horizontal  derivative  of  ^  is  retained  in  obtaining 
the  term  on  the  right  hand  side  of  (3.37). 

Before  deriving  a  parabolic  approximation  to  Eq.  (3.37),  the  need  for  such  an  approxi- 
mation should  be  discussed.  There  are  two  major  computational  drawbacks  to  numerically 
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solving  an  elliptic  equation.   The  first  is  that  solution  is  required  simultaneously  for  each 
grid.  This  means  that  for  a  large  computational  domain,  say  m  by  n  grids,  it  is  necessary 
to  invert  anmxnbymxn  matrix.   If  m  and  n  are  sufficiently  large,  problems  may  be 
encountered  in  terms  of  the  internal  memory  storage,  and  even  if  the  available  storage  is 
not  a  problem,  the  number  of  computations  required  to  invert  the  matrix  may  make  for 
an  exceedingly  slow  solution.  In  some  situations  this  may  be  acceptable;  however  for  the 
purposes  of  a  wave  induced  circulation  model  in  which  repeated  solution  for  the  wave  field 
is  necessary  the  program  would  take  a  very  long  time  to  run.  The  other  disadvantage  in- 
volved in  the  numerical  solution  of  an  elliptic  equation  is  that  boundary  conditions  must  be 
specified  at  all  of  the  boundaries.  In  many  applications  the  only  known  boundary  condition 
is  the  offshore,  or  incident  wave  amplitude.  For  example,  the  domain  of  interest  may  be  the 
wave  field  in  the  lee  of  an  obstruction  where  the  only  down  wave  condition  is  that  the  waves 
freely  pass  through  the  down  wave  boundary.  Or  in  the  present  study  in  which  the  wave 
amplitude  is  known  to  go  to  zero  at  the  shoreline  a  problem  ensues  in  that  the  shoreline  is 
(for  purposes  of  the  circulation  model)  at  the  grid  edge  and  not  at  the  grid  point. 

A  parabolic  equation  eliminates  both  of  these  problems.  Since  the  solution  proceeds  grid 
row  by  grid  row  where  each  solution  uses  the  results  of  the  previous  grid  row,  no  down  wave 
information  is  required.  Since  only  one  grid  row  is  solved  at  a  time,  the  solution  requires 
only  that  a  tridiagonal  matrix  equation  be  solved  to  obtain  values  for  the  grid  row.  The  only 
required  boundary  conditions  are  the  conditions  on  the  first  grid  row  (usually  the  offshore 
boundary)  and  lateral  boundary  conditions  which  are  usually  a  derivative  condition  that 
allows  for  the  wave  to  pass  through  the  boundary  without  any  reflection  at  the  boundary. 

There  are  several  ways  to  transform  the  elliptic  Eq.  (3.37)  to  a  parabolic  equation.  The 
most  direct  is  to  make  certain  assumptions  concerning  the  length  scale  of  variations  in  the  x 
direction.  At  the  heart  of  the  parabolic  approximation  is  the  assumption  that  the  direction 
of  wave  propagation  is  essentially  along  the  x  axis.  For  waves  propagating  at  an  angle  to 
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the  x  axis  the  proper  form  of  $  is 

I  =  -ig*i*±»leHfkco,td,+fk*n$dV)  (3  3 

and  the  proper  form  of  the  dispersion  relationship  is 

v  =  a  +  kcosOU  +  kainOV  (3.39) 

where  $  is  the  angle  of  the  wave  propagation  relative  to  the  x  axis.  The  essence  of  a  parabolic 
approximation  is  to  assume  that  the  angle  that  the  wave  rays  make  with  the  x  axis  is  small 
so  that  the  fcsinfl  term  can  be  neglected  in  Eqs.  (3.38)  and  (3.39)  and  cos0  is  assumed  to 
be  unity.  The  manner  of  doing  this  varies  with  different  authors.  Yan  (1987)  retains  the 
complete  dispersion  relation  (Eq.  (3.39))  but  deletes  the  y  dependency  in  Eq.  (3.38)  in  the 
derivation  of  his  governing  wave  equation.  He  solves  for  0  by  determining  the  slope  of  the 
instantaneous  water  surface.  For  long  crested  waves  this  method  is  satisfactory  although 
requiring  additional  computations.  However,  this  method  may  produce  erroneous  results  in 
diffraction  zones  where  the  presence  of  short  crested  waves  will  produce  a  horizontal  com- 
ponent of  the  slope  normal  vector  that  is  different  from  the  direction  of  wave  propagation. 
Kirby's  (1983)  method  is  to  make  the  assumption  that 

ju  Hfr^fej^**)  (340) 

so  that  the  etJk,lD'dv  part  of  the  phase  function  is  absorbed  into  the  amplitude  function, 
A. 

Using  Eq.  (3.40)  directly  in  Eq.  (3.37)  and  further  assuming  that  derivatives  of  Ax  are 
small  compared  to  derivatives  of  the  phase  function  (i.e.  that  ^  «  ikAx  )  will  yield  a 
parabolic  equation.  For  illustration  the  case  of  linear  diffraction  in  the  absence  of  currents 
and  in  a  region  of  constant  depth  is  examined.  Substitution  of  Eq.  (3.2)  into  the  Laplace 
equation  (3.6)  yields  the  Helmholtz  equation 

V2hi  +  k2i  =  0  (3.41) 
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This  equation  is  also  obtained  directly  from  Eq.  (3.37)  for  the  stated  assumptions  of  no 
currents,  no  dissipation,  and  constant  depth.    Substitution  of  (3.40)  into  the  Helmholtz 
equation  gives 

Axx  +  2ikAx  +  Avv  =  0  (3.42) 

and  employing  the  above  mentioned  assumption  concerning  derivatives  of  Ax  yields  the 
linear  parabolic  equation 

2ikAx  +  AyV  =  0  (3.43) 

Another  method  of  obtaining  a  parabolic  equation  is  to  use  a  splitting  method  in  which 
it  is  assumed  that  4>  is  composed  of  forward  moving  and  backward  moving  components  that 
are  uncoupled. 

$  =  i+  +  $-  (3.44) 

Uncoupled  equations  for  4>+  and  4>~  are  developed  and  since  the  equations  are  uncoupled 
only  the  equation  for  4>+  is  retained.  In  the  following  the  work  of  Kirby  (1983)  and  Booij 
(1981)  is  closely  followed. 

The  second  order  elliptic  equation 


can  be  split  exactly  into  equations  governing  forward  and  backward  disturbances  $+  an  $~ 
which  are  given  by 

"aT  "  +t^  (3-46) 


Booij  chose  the  relation 


(3.48) 


*  =  W  (3.49) 
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where  £  is  an  as  yet  unknown  operator  on  $.  Substituting  into  Eq.  (3.45)  and  neglecting 
derivatives  of  £  alone  leads  to  the  result 

*~+(f)       (^Jk  +  -r2$  =  0  (3.50) 

The  elliptic  equation  (3.37)  is  now  put  into  the  form  of  Eq.  (3.50).  Letting  p  =  CCg 
for  notational  convenience  and  expanding  some  of  the  terms  of  Eq.  (3.37)  in  their  spatial 
coordinates  gives 

{p<f>z)z  +  {p4>v)v  +  k2pj>  +  (u/2  -  a2  +  i'w(Vfc  •  U))4>  +  iaw4> 
~{U2k)z  ~  (V24>v)v  -  [UV$t)v  -  [UVfai  +  2iuU  -Vj  =  0  (3.51) 

Booij  chose  to  neglect  terms  involving  the  square  of  the  mean  current  components  assuming 
that  they  are  smaller  than  terms  involving  p  and  then  arranged  the  above  equation  as 

kz  +  P_1(p«  +  2ivU)$t  +  *2  (l  +  pr)  I  =  0  (3.52) 

where 

Ml  =  (w2  -  a2  +  iu(Vh  •  U)  +  i<rw)i  +  2iu;V4v  +  (p^,)„  (3.53) 

Kirby  points  out  that  this  choice  by  Booij  leads  to  his  inability  to  derive  the  correct  wave 
action  equation.  Kirby  instead  does  not  drop  any  of  the  squares  of  the  mean  current 
components  until  after  the  splitting  is  complete.  By  assuming  that  the  waves  are  oriented 
in  the  x  direction  the  dispersion  relation  is  written  as 

a  =  v  -  kU  (3.54) 

leading  to  the  result 


2  2 


-  a2  =  2wkU  -  (kU)2  (3.55) 

Using  this  relationship  Kirby  then  arranged  Eq.  (3.51)  as 

kz  +  (p  -  V*)-l(p  -  U2)Jt  +  k2  (l  +  k2{p^u2))  l  =  0  (3.56) 
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where 


+  {(P-^2K}V  +  2.W.V^  (s.57) 

Comparison  with  (3.50)  gives 

i1  +  k*(p"u*))  <3-58) 


72    =    **(!  + 


e2 


-      =    (P-U)  (3.59) 

or  solving  for  7  and  £ 

Booij  points  out  that  the  splitting  matrix  approach  employed  by  Radder  (1979)  is  equivalent 
to  the  choice 

Using  the  results  (3.60)  and  (3.61)  in  Eq.  (3.46)  along  with  the  definition  (3.49)  gives  the 
parabolic  equation 

By  expanding  the  pseudo-operators  in  the  form 

(1+t^j)^+     -     (I  +  iP(^j)*+  (363) 

the  result  of  both  expansions  may  be  summarized  in  the  general  equation 
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where 

Radder 


»-! 


Bboij  (366) 

and  where  Pi  =  P\  +  \  for  both  approximations.  Kirby  chose  the  approximation  with 
Pi  =  0  since  for  his  derivation  which  was  for  weakly  nonlinear  waves  it  should  correspond 
to  the  information  contained  in  the  nonlinear  Schrodinger  equation.  With  Pi  =  0  Eq.  (3.65) 
becomes 

*(p  "  U2)fc  +  1  {k{p  -  U2)}z  p  =  ik\p  -  U7)4>+  +  %-M4>+  (3.67) 

where  M^+  is  defined  by  Eq.  (3.57).  An  equation  for  the  complex  amplitude  A  is  then 
obtained  by  making  the  substitution 

j>+  =  -ig(^)eifkdx  (3.68) 

which  after  dropping  squares  of  the  components  of  the  mean  current  results  in 

Kirby  then  further  modified  the  equation  by  a  shift  to  a  reference  phase  function  using  the 
relation 

A  =  A'ei(ix-fkd')  (3.70) 

where  k  can  be  taken  as  an  averaged  wave  number,  resulting  in 

{C,  +  U)A'X  +  .•(*  -  k)(C,  +  U)A'  +  °-  (^^)    A' 


•  *    L  X  "    '  VI  y 


The  shift  to  a  reference  phase  function  is  intended  to  incorporate  the  phase  information 
in  the  complex  amplitude  so  that  the  instantaneous  water  surface  can  be  easily  recovered 
in  the  form 

t7(x,y)  =  »{AV/idl}  (3.72) 
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where  »  designates  the  real  part.  This  is  an  important  consideration  especially  when  ob- 
taining the  radiation  stresses. 

It  is  important  that  in  obtaining  Eq.  (3.69)  that  squares  of  the  products  of  the  current 
components  be  retained  until  the  last  step.  The  reason  for  this  is  seen  in  that  terms  such 
as  WV  and  kU*  are  combined  with  terms  of  uV  and  uU  using  the  dispersion  relationship 
to  give  terms  of  aU  and  aV . 

Equation  (3.71)  has  been  successfully  used  by  Kirby  in  the  form  given  and  with  modifi- 
cations. The  modifications  are  the  addition  of  a  non-linear  term  which  makes  the  equation 
applicable  for  weakly  non-linear  waves.  However  for  present  purposes  Eq.  (3.71)  presents  a 
certain  problem  which  will  shortly  be  corrected.  The  problem  is  that  for  waves  propagating 
at  an  angle  to  the  x  axis  it  is  not  possible  to  obtain  a  correct  form  of  the  conservation 
of  wave  action.  The  case  of  no  current  on  a  beach  with  straight  and  parallel  contours  is 
given  as  an  illustration.  The  equation  for  the  conservation  of  wave  action  reduces  to  the 
conservation  of  wave  energy  flux.  In  mathematical  terms  this  is 

[CgCos0a2]^  =  O  (3.73) 

where  a  is  the  wave  amplitude.  For  waves  approaching  the  beach  in  a  direction  normal  to 
the  beach  Eq.  (3.71)  reduces  to 

C,<+\(C„)xA'  =  0  (3.74) 

Multiplying  Eq.  (3.74)  by  the  complex  conjugate  of  A',  and  add  to  it  the  complex  conjugate 
equation  of  (3.74)  multiplied  by  A',  yields 

fcl<]x  =  0  (3.75) 

which  is  Eq.  (3.73)  for  the  special  case  of  cos*  =  1.  However  it  is  not  possible  to  obtain 
Eq.  (3.73)  from  Eq.  (3.71)  for  the  general  case  of  waves  at  angle  6  on  a  beach  with  straight 
and  parallel  contours. 

The  importance  of  this  inability  to  give  the  correct  form  of  Eq.  (3.73)  is  that  if  energy 
flux  outside  the  breaker  line  is  not  conserved  then  gradients  of  the  Sxy  radiation  stress 
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will  generate  currents  off-shore  of  the  breaker  line  in  the  circulation  model.  Although 
such  currents  in  general  will  be  of  a  small  magnitude,  they  can  cause  some  numerical 
computational  problems.  The  problem  though  is  easily  corrected  by  using  as  the  dispersion 
relationship 

w  =  a  +  k  cos  9  U  (3.76) 

and  the  substitution 

^=_.^y)e,(/tCo.^x)  (3?7) 

Rather  than  using  a  splitting  method,  the  procedure  will  be  to  go  directly  from  Eq.  (3.51) 
using  Eq.  (3.77)  and  writing 

w2  -a1  =  2ukcos9U  -(kcos9U)*  (3J8) 

in  place  of  Eq.  (3.55).  Dropping  the  (£)«  term  and  all  terms  containing  products  of  the 
current  components  that  cannot  be  reduced  as  explained  above  results  in  the  following 
parabolic  equation. 


(CgcoB9  +  U)Ax  +  - 


Cg  cos  9  +  U 


A  +  VAv+-(-\    A 


Jx  -         /f 


m. 


-\kCg{l-cos29)A-i-\cCg{^    +^A    =    0  (3.79) 

Now  adding  a  phase  shift 

A  =  A,eiUicotSd'-fkco"d*)  (3  80) 

so  that  the  free  surface  can  be  recovered  by 

r,(x,y)  =  A'eiUl™id*)  (3.81) 

where  9  is  the  wave  angle  obtained  by  Snell's  law  in  the  absence  of  any  obstructions  or  y 
variation  of  the  depth  or  current.  The  following  governing  equation  is  thus  obtained. 

(C,  cos  9  +  U)A'X  +  i(k  cos  9  -  k  cos  9){C9  cos  9  +  U)A'  +  -  [C'cosg  +  tH    A> 

+VA'< +5(7)/-  s*wi  - cos2  w  -  i  \cc°(tA  +  V = °  (3-82^ 


45 
Details  of  the  algebra  and  assumptions  involved  in  obtaining  Eq.  (3.79)  are  given  in  Ap- 
pendix D.  For  computational  purposes  it  will  be  assumed  that  cos  9  =  cos  9. 

From  Eq.  (3.82)  the  correct  form  of  the  conservation  of  wave  action  can  be  obtained. 
For  the  case  of  no  y  variation  of  the  wave  number  Eq.  (3.82)  reduces  to 

C,  cos0  +  W 


(C,*»9  +  U)A'm  +  Z 


\*+**+m* 


-i*C,(l-cos2*)A'-i[cC,(^),,]    +  |a    =    0         (3.83) 

Multiplying  Eq.  (3.83)  by  the  complex  conjugate  of  A  and  adding  to  it  A  times  the  complex 
conjugate  equation  of  (3.83)  yields 


{Ctcos9  +  U) 


+ 


■Ml 


=  o 


(3.84) 


which  is  a  closer  approximation  to  the  exact  expression  which  should  have  a  (Ct  sin  9  +  V) 
term  instead  of  V  in  the  y  derivative  term. 

33     Energy  Dissipation  Due  to  Wave  Breaking 

Kirby  (1983)  found  a  method  of  relating  the  dissipation  coefficient  w  to  the  work  of 
Dally  (1980)  who  proposed  that  the  decay  of  energy  flux  in  the  surf  zone  was  proportional 
to  the  amount  of  excess  energy  flux.  The  excess  energy  flux  is  the  difference  between  the 
actual  flux  and  a  stable  energy  flux.  In  mathematical  terms  this  relationship  is  expressed 

as 


K 


(EC9)z  =  -TlEC>-(ECo)'l 


(3.85) 


where 


ECg    =    the  energy  flux 

E    =    wave  energy  =  -pgH2 

8 

Ct    =    the  rate  of  energy  propagation  which  is  — 

dk 

K    —    constant 


h    =    the  water  depth 
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(ECg),    =    the  'stable'  energy  flux. 

H    =    the  local  wave  height 

The  stable  energy  flux  is  obtained  by  assuming  it  to  be  a  function  of  the  height  of  a  wave 
propagating  over  a  flat  bottom  after  breaking  ceases,  and  that  this  height  is  a  constant 
times  the  water  depth,  (f  )^  =  7.  The  constants  K  and  7  were  determined  by  analysis  of 
experimental  data  obtained  from  Horikawa  and  Kuo  (1966). 

Kirby  related  Daily's  work  to  the  dissipation  term  w  by  deriving  a  conservation  equation 
for  the  wave  action  f .  Garrett  (1967)  and  Bretherton  and  Garrett  (1968)  showed  that  for 
waves  propagating  in  a  moving  and  variable  medium,  the  quantity  that  is  conserved  is  the 
wave  action  (f  J ,  which  is  the  wave  energy  divided  by  the  intrinsic  frequency.  Kirby  derived 
a  wave  action  conservation  equation  from  the  linear  hyperbolic  wave  Eq.  (3.34)  using  the 
definition  of  W  as  given  by  Eq.  (3.35): 

D26  -  Dd> 

^JT+  (V*  •  &)jZ  -  Vk(CC,Vk4)  +  (<rJ  -  k2CCt)$  =  iW*  (3.86) 

Kirby  follows  the  approach  of  Jonsson  (1981)  assuming  that  4>  can  be  represented  as 

i>  =  Reie  (3.87) 

where  R  and  6  are  real.  From  linear  theory  it  is  obvious  that 


R=^-\A\ 


a 


To  lowest  order  the  derivatives  of  6  are  approximated  as 

6«  =  -u 


(3.88) 


(3.89) 


Vfc©  =  *  (3.90) 

Substituting  into  (3.86)  gives  the  complex  equation 

D2R      ,„      ~,DR 

-^2  +  (V*  •  U)-jh  ~  Vh{CC„VhR)  -  i(atR  +  2aRt)- 

«'  { V*  •  (Ua)R  +  2{aU)  •  VhR  +  V„  •  (kCC9)R  +  2kCC9  •  Vkr) 

-iawj>    =    0        (3.91) 
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Since  R  and  6  are  real,  the  real  and  imaginary  parts  of  (3.91)  must  individually  equal  zero. 
Taking  the  imaginary  part  and  multiplying  by  R  gives 

{oR2)t  -  *0  •  VhR2  -  V*  •  (aff)R2  -  CC9k  ■  VkR2 

-Vk.{kCCt)R*-<rwR2    =    0  (3.92) 

which  using  the  substitution  kCCg  =  aCg  is  written  as 

{aR2)t  +  V„  •  (oR2{U  +  C9j)  =  -awR2  (3.93) 

Using  (3.88)  R2  may  be  replaced  by 

which  yields  the  wave  action  conservation  equation 

(?).+M(?)«»+*}"(5)  m 

Assuming  a  time  steady  wave  field  and  neglecting  currents,  (3.95)  becomes 

(£C,)2  =  -wE  (3.96) 

Comparing  (3.96)  with  (3.85)  and  noting  that  (Ct).  =  C„  w  is  written  as 


w=^^ 


where  ^T  is  the  local  wave  height  given  by  2\A\. 

In  obtaining  Eq.  (3.97)  Kirby  neglected  the  effects  of  the  currents.  However,  the  same 
expression  (3.97)  is  obtained  when  currents  are  included  in  the  derivation.  Dally  (1987) 
extended  his  work  to  investigate  wave  breaking  on  currents  and  proposed  the  following 
equation: 

[(C*  +  V){§)l  =  ~[EC,-(EC,).)  (3.98) 

For  the  case  of  a  steady  state  wave  field  propagating  in  the  x  direction  the  conservation  of 
wave  action  Eq.  (3.95)  reduces  to 
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Setting  the  right  hand  sides  of  Eqs.  (3.98)  and  (3.99)  equal  to  each  other  results  in  Eq.  (3.97). 
At  first  glance  it  may  appear  unusual  that  the  same  expression  for  the  dissipation 
coefficient  results  for  the  case  of  currents  as  for  the  case  of  no  currents.  This  can  be 
explained  in  part  by  the  fact  that  Daily's  work  was  with  waves  where  breaking  was  depth 
limited.  For  the  case  of  current  limited  breaking  (for  example  a  stopping  current),  another 
expression  is  necessary.  In  this  study  Eq.  (3.97)  will  be  used,  since  in  general  the  current 
field  will  be  moderately  weak.  Such  an  assumption  of  weak  currents  was  invoked  in  the 
derivation  of  the  governing  parabolic  equation. 


CHAPTER  4 
FINITE  DIFFERENCE  SCHEMES  AND  BOUNDARY  CONDITIONS 


4-1  Numerical  Solution  of  the  Governing  Equations 
The  finite  differencing  of  the  governing  equations  in  the  circulation  model  is  derived 
by  a  matrix  analysis  as  suggested  by  Sheng  and  Butler  (1982).  A  rational  analysis  of  the 
complete  Eqs.  (2.1-2.3)  including  all  the  nonlinear  terms  is  not  possible.  However  guidance 
can  be  obtained  through  an  analysis  of  the  linearized  equations.  For  this  purpose  the  linear 
equations  for  the  special  case  of  a  flat  bottom  and  without  any  of  the  forcing  terms  will  be 
examined. 

dfj       _dU       _dV 

(4.1) 

(4.2) 
Let  W  be  a  vector  of  fj,U  and  V 


dfj      au    Ddy__ 

dt  dx  dy 

dU        dfj 

dV        dfj 


w  = 


n 

u 

V 


(4.4) 


Using  Eq.  (4.4),  the  governing  equations  are  written  as 


Wt  +  AWt  +  BWv  =  0 


(4.5) 


where  A  and  B  are  the  following  matrices 


A  = 


0 

D    0  ' 

9 

0     0 

0 

0     0 

and         B  = 


0 

0 

D  ' 

0 

0 

0 

9 

0 

0 

(4.6) 


To  find  an  implicit  finite  difference  scheme  for  the  governing  equations,  Eq.  (4.5)  is 
finite  differenced  as 

Wn+1  -  Wn 

— +  AW^  +  BW?+1  =  0  (4.7) 
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This  equation  is  solved  for  Wn+1  in  terms  of  Wn 

(1  +  2XX  +  2A„)Wn+1  =  Wn  (4.8) 

where 

1     AT  1     AT 

A*=2AAX**         and        A"=2SAY^  (4'9) 

AXxXvWn+1  is  added  to  the  left  hand  side  of  the  equation  and  2XyWn  is  added  and  sub- 
tracted to  the  right  hand  side  of  Eq.  (4.8)  so  that  the  equation  can  be  factored  to  get 

(1  +  2A,)(1  +  2Xy)Wn+1  =  (1  -  2Xv)Wn  +  2XyWn  (4.10) 

An  intermediate  value  W*  is  introduced  so  that  the  above  equation  can  be  written  as 
two  equations 

{l  +  2Xx)W*  =  {l-2Xy)Wn  (4.11) 

(1  +  2Xv)Wn+1  =  W*  +  2XvWn  (4.12) 

Note  that  by  eliminating  W*  from  the  two  Eqs.  (4.11)  and  (4.12),  Eq.  (4.8)  is  recovered 
with  the  addition  of  an  error  of  AXxXyWn+l  which  has  been  introduced. 
Writing  out  the  component  equations  of  Eq.  (4.11)  gives 

<t*  +  *>^W  =  tl«  -  D^6yV»  (4.13) 

U*+g^8xfj*  =  Un  (4.14) 

Vt=Vn-9^yr  (4.15) 


AY 


(4.16) 


While  writing  out  the  component  equations  of  Eq.  (4.12)  gives 


AT  AT 

/n+1  +  D—SyV^1  =  fj*  +  D^6VV»  (4.17) 


Un+1  =  U*  (4.18) 

AT 

AY""''  '      rffAY' 


vn+1  +  ^VT1  =  V  +  g^W  (4.19) 
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Using  Eq.  (4.15)  to  eliminate  V*  in  Eq.  (4.19)  and  (4.18)  to  eliminate  U*  in  Eq.  (4.14)  the 
following  are  obtained.  For  the  x  direction 

AT  AT 

^  +  DAXS*U"+1  =  *"  "  D^Vn  (4-20) 

AT 
Un+1  +  g~Sxfj*  =  Un  (4.21) 

and  for  the  y  direction 

**1  +  D^n+1  =  V*  +  D^6yV»  (4.22) 

AT 
^ +  #££*f"  =  V  (4.23) 

The  above  equations  are  numerically  stable  for  any  combination  of  AT,  AX,  and  AY. 
However  as  a  circulation  model  they  are  incomplete  in  that  they  do  not  contain  any  of  the 
forcing  terms  that  drive  the  circulation.  In  particular  in  this  model  wave  forcing  represented 
by  gradients  in  the  radiation  stresses  and  bottom  friction  are  of  concern.  Bottom  friction  is 
especially  important  because  at  steady  state  conditions  gradients  in  the  radiation  stresses 
must  be  balanced  by  either  gradients  in  the  water  surface  elevation  or  by  bottom  friction 
forces.    In  the  absence  of  friction  forces  the  model  could  yield  continuously  accelerating 
currents.  The  method  of  adding  the  forcing  terms  is  somewhat  arbitrary  in  that  the  rational 
derivation  of  the  finite  differencing  scheme  of  Sheng  and  Butler  cannot  be  expanded  to 
include  these  essentially  non-  linear  terms.   The  addition  of  these  terms  to  the  model  also 
limits  the  latitude  of  the  ratio  of  |£  in  a  manner  that  is  not  analyzable.   That  both  the 
friction  and  the  radiation  stresses  are  nonlinear  should  be  recognized.    The  non-linearity 
of  the  velocity  friction  term  is  quite  evident  since  it  is  represented  by  the  quadratic  term 
pF\U\U.    The  radiation  stress  terms,  which  represent  excess  momentum  flux  induced  by 
waves,  are  functions  of  the  wave  height  squared,  based  on  linear  wave  theory.   In  shallow 
water  the  wave  height  H  is  a  function  of  the  total  water  depth,  h0  +  fj  so  that  H2  is  a 
function  of  fj2. 

All  the  forcing  terms  and  the  nonlinear  acceleration  terms  are  added  at  the  present  or 
n*    time  level  with  the  exception  of  the  bottom  friction  term.  The  quadratic  friction  term 
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is  finite  differenced  as  pF\Un\Un+1  and  the  friction  term  as  suggested  by  Longuet-Higgins  is 
differenced  as  pF  \u„h\  Un+1,  where  |uort|  is  in  terms  of  values  at  the  n  time  level.  However, 
when  using  the  integral  quadratic  expression,  it  is  evaluated  entirely  in  terms  of  the  present 
n      or  known  values. 

The  nonlinear  advective  acceleration  terms  are  added  at  the  nth  time  level  using  central 
difference  derivatives.  In  the  x  direction  momentum  equation  these  terms  are 

U'^'iAX~"  +  i«#  +  *■«  +  HW  +  V<-W)^y^-'  (4.24) 

and  in  the  y  direction  equation  they  are 

\{utJ  +  u^j  +  cr<i+1  +  tWO1**1^-1'  +  K^Sfli^td         (4.25) 

In  the  model  of  Vemulakonda  (1984)  which  is  based  upon  the  WIFM  (Waterways  Ex- 
periment Station  Implicit  Flooding  Model)  (Butler  1978)  model,  a  three  time  level  leap-frog 
scheme  is  used.  This  gives  the  advantage  that  all  of  the  forcing  terms  (and  the  nonlinear 
acceleration  terms)  are  added  at  the  time  centered  level.  However  this  has  some  disad- 
vantages in  that  additional  computer  memory  and  computational  time  are  required  and 
leap-frog  schemes  can  exhibit  instabilities.  In  the  present  study  a  two  time  level  scheme 
is  used  and  the  wave  forcing  terms  as  well  as  the  non-linear  acceleration  terms  are  added 
at  the  present  time  level.  Since,  for  all  of  the  applications  in  this  report,  the  model  is  run 
until  steady  state  conditions  are  obtained,  it  is  reasoned  that  there  is  no  great  necessity  to 
properly  time  center  all  of  the  terms.  If  the  use  of  the  model  is  extended  to  time  dependent 
applications,  this  may  need  to  be  corrected. 

4.2     Finite  Differencing  of  the  Parabolic  Wave  Equation 
Dividing  Eq.  (3.82)  by  (C,cos0  +  U)  yields 

'  2(C,cos0  +  tOA  +  (C,cos*  +  CO    * 

+ z (L)  a'    IBSzggflv 

2(Ctcosff  +  U)  \a)y  2  (C,cos*  +  tf)  A 

-2(C,cos*  +  tf)  [CC'^)V  +  2(Ctcos9  +  U)A'  =  °  (4"M) 
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The  following  notation  is  used  to  describe  the  finite  differencing  of  Eq.  (4.26):  F*- 
represents  the  value  of  F(i,j)  where  F  is  any  term  and  i  and  j  designate  the  grid  row 
and  grid  column.  Because  of  the  staggered  grid  system  used  for  the  circulation  model  as 
previously  explained,  in  which  the  velocity  components  are  specified  at  the  grid  edges,  the 
velocity  terms  U  and  V  are  averaged  across  the  grid  so  as  to  produce  a  grid  centered  term 
for  use  in  the  wave  model.  Thus  Uf  in  the  wave  model  is  actually  \(Ui.  +  Ui+l)  »nd  similarly 
for  V  and  other  grid  designations. 

A  Crank-Nicholson  finite  differencing  scheme  is  used  for  Eq.  (4.26).  All  terms  not 
involving  an  x  derivative  are  an  average  of  the  value  of  the  term  on  the  i  grid  row  and  on 
the  i  +  l  grid  row.  Terms  having  y  derivatives  are  centered  at  i ,  j  on  the  i  row  and  i  + 1 ,  j 
on  the  i  +  l  grid  row.  The  one  term  with  an  x  derivative  is  written  as  ^y^1 .  Thus  all 

AA 

terms  are  centered  at  i  +  },/.  The  finite  differencing  of  Eq.  (4.26)  is  written  as 

AX 

+  l-  [(*'+1cos*'+1  -  ^eo.^1)^1  +  (jfc'cos**  -  *j.cos?')A«] 

■  i  i  w« cos °  ±  u)H+l  -  i(gp  «* ; + u)/°}\  i  / ,- . ,  aA 

2 AX  \  [(Ct  cos  6  +  U)/o])+i  +  {(Cg  cos  0  +  U)/a})  j  ^     +  Ai) 


1 

+  2 


Va 


(Ca  cos  0  +  U) 


).  2AY       -+[ 


(Cgcos9  +  U) 


2AY 


+  2  1  \2(CtcoB0  +  U)l  . 
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( i__V 

\2(CtcoB0  +  U)J  . 


2AY 
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*i 
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fcCg(l-cos2fl) 
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JfcC,(l-cos20) 


■(*) 
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■    I+l 

v. 
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(*). 
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i      \ 


4    (Cg  cos  e  +  u)*1  T  (Cg  cos  e  +  u)y 
1/     w^a^1  w^a;       ) 

4  \  (Cgcos0  +  r/)j+»  +  (Cgcosfi  +  U))  f~°    (427) 
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where  p  =  CCt  for  notational  convenience  and  the  y  derivative  of  p  (£)  is  written  as 

In  order  to  evaluate  wj+J  by  the  method  to  be  explained  in  section  3.3  an  intermediate 
value  of  Af1  is  obtained  by  the  explicit  formulation 


+ 


( 


_^ Vf 

2{Cgcos9  +  U)J  .  [  2AY 


A'"1"1  -  A' 
'  AX    '  +  i{V  cos  0'  -  k)  cos  0')AJ 

[(C,  cos  9  +  V)/o^1  -  [{Cg  cos  9  +  f/)/g]»  )      . 
AX  I  [(C,  cos  0  +  CT)/flrj}+l  +  [(C,  cos  0  +  r/)/<r]j  J  A%i 
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2  (C,  cos  9  +  tf)$  +  2(C,  cos  $'+ U))  =°      (429) 
When  Eq.    (4.27)  is  written  out  using  Eq.  (4.28)  and  w*+1  obtained  through  use  of 
(4.29)  there  results  for  any  particular  values  of ."  and  j  an  equation  for  the  three  unknown 
values  ilj+l,  A)+1,  and  A)t\.  The  equation  is  put  into  the  form 

aiA^\  +  bjA^  +  cjA^\  =  dj  (4.30) 

where  the  coefficients  Qj,  bjt  and  Cj  are  known  quantities  involving  the  properties  of  the 
environment  (depth,  current,  celerity,  etc.)  and  the  dj  involves  the  environmental  properties 
and  the  known  values  of  the  amplitude  on  the  i  grid  row.  When  Eq.  (4.30)  is  written  for 
all  the  j  values,  j  equal  to  1  through  N,  for  a  specified  ,"  value,  the  ensemble  of  N  equations 
takes  the  form  of  a  tridiagonal  matrix  equation.  For  illustrative  purposes  let  N  be  equal  to 
four.  The  resulting  matrix  equation  is 


a\     6i     ci 
ai     62     C2 
<*s    h    cs 

<*4      *4      C< 


r4,+i- 

A\+1 

■  <*i  ] 

A+1 

d2 

A+1 

d3 

4+1 

d4 

.  *s+1  J 

(4.31) 
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Equation  (4.31)  contains  four  simultaneous  equations  with  six  unknowns.  Two  bound- 
ary conditions  should  be  specified.  At  each  boundary  a  functional  relationship  is  established 
between  the  grid  inside  the  domain  and  the  grid  outside  of  the  domain.  For  the  wave  prop- 
agation model  there  are  two  types  of  boundary  conditions  that  are  readily  available  for  use 
in  a  computer  model.  These  are  a  totally  reflective  boundary  condition  which  would  be 
used  to  model  a  laboratory  experiment  in  a  wave  basin,  and  non-reflecting  non-interfering 
boundary  which  will  allow  waves  to  freely  pass  through  the  boundary.  For  the  reflective 
boundary,  the  mathematical  statement  of  the  boundary  condition  is  |^  =  0.  At  the  lower 
boundary  this  is  finite  differenced  as  ^^  =  0  which  gives  the  functional  relationship 

M  =  M-  (4.32) 

At  the  upper  boundary  the  functional  relationship  for  a  reflective  boundary  is 

AN+i  =  AN.  (4.33) 

For  the  non-reflective  boundary  the  mathematical  statement  is  |^  =  -imA  where  m  is 
the  longshore  component  of  the  wavenumber.  Since  the  lateral  boundaries  are  sufficiently 
far  from  any  diffraction  effects,  only  refraction  and  shoaling  effects  are  present.  Thus 
using  Snell's  Law  m  =  ksmO  =  kosm0o.  The  non-  reflective  boundary  condition  is  finite 
differenced  at  the  lower  boundary  as 

A\  -  A0       im 

AY      =  ~Y  lAi  +  M  (4.34) 

which  gives 

Similarly  at  the  upper  boundary  the  functional  relationship  used  to  eliminate  AN+1  is 

AN+1  =    h        imfY\A»  (4"36) 

The  use  of  either  pair  of  boundary  conditions,  (4.32)  and  (4.33)  or  (4.35)  and  (4.36)  to 
eliminate  A0  and  AN+1  will  then  yield  an  N  by  N  tridiagonal  matrix  equation  comprising  N 
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Figure  4.1:  Definition  sketch  of  grid  nomenclature  and  coordinate  system. 

equations  for  N  unknowns.  Returning  to  the  illustrative  example  where  N  is  equal  to  four, 
the  matrix  Eq.  (4.31)  becomes 


h  +  7\&\     ci 

fl2  ^2       ^2 

a*    bt  +  Tnca 


\  Al+1 1 

\di] 

A\+i 

d2 

4+1 

dz 

A\+1 

d4 

(4.37) 


where  ?\  and  Tn  represent  the  functional  relationships  as  described  above  on  the  j=\  and 
j=N  boundaries. 

Equation  (4.37)  is  solved  easily  for  the  unknown  Aj+1  using  a  double  sweep  method. 
In  the  computer  program  the  subroutine  CTRIDA  uses  an  algorithm  of  Carnahan,  Luther, 
and  Wilkes  (1969)  to  efficiently  accomplish  the  double  sweep  solution  of  eq.  (4.37). 
4.3     Boundary  Conditions  in  the  Circulation  Model 

As  shown  in  figure  4.3  a  staggered  grid  system  is  used  so  that  velocities  are  expressed 
at  the  grid  edge  and  fj  and  the  wave  properties  are  defined  at  the  grid  center.  So  that  the 
wave  model  will  have  full  grids  it  is  desirable  to  place  the  boundaries  of  the  model  at  the 
grid  edges  on  all  sides,  with  the  exception  of  the  offshore  grid  row  where  the  wave  conditions 
are  an  input  or  boundary  condition  for  the  wave  model.  Figure  1.1  can  also  be  referred  to. 

The  offshore  boundary  condition  that  is  used  in  the  circulation  model  is  that  the  mean 
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water  surface  level  is  not  changed  by  any  of  the  circulation  within  the  model,  mathematically 

ETA(1,J)  =  0  J  =  i,N  (43g) 

This  condition  necessitates  that  the  offshore  boundary  be  sufficiently  far  from  any  set-up  or 
set-down.  Other  possible  offshore  conditions  are  a  no-flow  condition  or  a  radiation  condition. 
The  no  flow  condition  would  conserve  water  within  the  computational  domain  so  that  the 
volume  of  water  involved  in  set-  up  would  be  compensated  by  an  equal  volume  of  set-down 
offshore.  In  an  open  ocean  environment  the  open  ocean  can  be  conceived  of  as  an  infinite 
reservoir.  However  in  modeling  laboratory  experiments  a  no  flow  condition  as  the  offshore 
boundary  condition  would  be  the  most  appropriate.  A  radiation  boundary  condition  as 
proposed  by  Butler  and  Sheng  (1984)  imposes  the  condition 

dfi       „  dfj 

*+C£  =  °  (4-39) 

This  condition  allows  long  waves  to  freely  propagate  out  of  the  computational  domain 
without  any  reflection.  This  tends  to  eliminate  any  seiche  within  the  computational  domain 
which  is  induced  by  the  start-up  transients.  This  condition  was  applied  at  one  time  in  the 
present  model  but  it  did  not  preserve  the  offshore  zero  set-up  condition  and  was  changed  in 
favor  of  the  more  rigid  zero  water  surface  elevation  condition.  If  the  model  is  adapted  for 
use  with  time  varying  wave  trains  (an  adaptation  that  the  model  can  easily  accommodate) 
then  the  radiation  condition  will  become  imperative. 

At  the  shoreline  the  boundary  condition  is  that  there  is  no  flow  into  or  out  of  the  beach. 
Since  the  model  allows  for  flooding  of  the  beach  the  index  of  the  shoreline  changes.  The 
model  allows  for  differential  flooding  along  the  shoreline,  for  example,  due  to  differences  of 
wave  intensity  in  the  lee  of  structures.  For  each  grid  column  in  the  y  direction  the  last  grid 
is  the  IWET(J)  grid,  where  IWET(J)  is  allowed  to  change  with  changes  of  the  mean  water 
level.  The  numerical  statement  of  the  no  flow  boundary  condition  at  the  shoreline  is 

C/(IWET(J)+1,J)  =  0  j=i)N  (4  40) 
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This  boundary  condition  is  used  in  other  models,  though  not  all  models  allow  for  differential 
flooding  of  the  shoreline.  Yan  (1987),  whose  model  was  used  to  investigate  waves  and 
currents  around  a  tidal  or  fluvial  discharge,  used  the  no  flow  condition  away  from  the 
discharge  and  defined  tf(IWET(J)+l,J)  as  equal  to  the  discharge  per  grid  divided  by  the 
water  depth  at  the  grids  where  discharge  occurred. 

Other  conditions  are  also  imposed  at  the  shoreline  boundary.  These  are  that  the  wave 
height  and  the  longshore  velocity  are  both  zero.  These  conditions  are  imposed  at  the 
shoreward  grid  edge  of  the  IWET(J)  grids  even  though  the  wave  height  and  longshore 
velocities  are  specified  a  half  grid  space  form  the  shoreline.  This  is  done  by  writing  x 
derivatives  of  the  radiation  stress  terms  and  of  the  V  velocities  that  need  to  be  positioned 
at  the  center  of  the  IWET(J)  grid  in  terms  of  three  points.  The  three  points  are  at  the 
IWET(J)-1  and  IWET(J)  grids  and  the  zero  value  at  the  shoreline  at  the  edge  of  the 
IWET(J)  grid. 

For  the  lateral  boundary  conditions  several  alternatives  are  available.  Birkemeier  and 
Dalrymple  (1976)  and  Ebersole  and  Dalrymple  (1979)  both  used  a  periodic  or  repeating 
beach.  This  meant  that  computations  were  performed  between  the  second  and  Nth  grid, 
and  values  outside  of  the  computational  domain  are  determined  according  to 

Q(I,1)  =  Q(I,N)  (441) 

q(i,n  +  i)  =  q(i,2)  foraln  (442) 

for  any  quantity  Q. 

At  one  point  this  was  used  in  the  present  model  until  the  problem  of  what  should  have 
been  transient  edge  waves  generated  by  wave  obstructions  during  the  start-up  were  found 
to  be  trapped  in  the  model.  The  periodic  boundary  condition  allowed  for  disturbances 
propagating  out  of  the  domain  at  one  lateral  boundary  to  enter  the  domain  at  the  other 
boundary.  For  straight  and  parallel  contours,  the  periodic  beach  condition  worked  well  in 
the  present  model  but  with  the  presence  of  breakwaters  another  condition  had  to  be  used. 
For  the  lateral  boundary  conditions  there  are  two  choices,  either  closed  or  open  bound- 
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aries.  Closed  boundaries  are  easily  applied  by  the  condition  of  no  flow  across  the  boundary. 
Such  a  choice  of  boundary  conditions  would  accurately  duplicate  the  physical  conditions  ex- 
isting in  a  laboratory  basin  (or  for  prototype  conditions  within  a  harbor).  Open  boundaries 
are  obtained  by  specifying  derivative  conditions.  In  the  present  study  an  open  boundary 
is  obtained  by  specifying  that  there  are  no  gradients  of  the  depth  and  the  wave  proper- 
ties across  the  lateral  boundaries.  This  necessitates  that  the  lateral  boundaries  are  placed 
sufficiently  far  from  any  obstructions  or  any  y  gradients  caused  by  the  obstructions.  An 
alternative  open  boundary  condition  would  be  that  there  is  no  gradient  of  the  flow.  This 
would,  however,  put  the  position  of  the  boundary  condition  at  the  center  of  the  bordering 
grid  rather  than  at  the  edge  of  the  grid  and  could  possibly  introduce  error  at  the  boundary. 
4.4     Boundary  Conditions  in  the  Wave  Model 

In  the  wave  model  only  lateral  boundary  conditions  are  needed  for  computational  pur- 
poses. The  offshore  condition  is  merely  the  specified  wave  amplitude  and  wave  angle.  No 
condition  is  needed  at  the  shoreline  since  the  parabolic  wave  model  marches  grid  row  by 
grid  grow  and  needs  only  the  lateral  conditions. 

The  offshore  condition  is  that  each  time  the  wave  model  is  called  the  wave  amplitude 
on  the  first  grid  row  is  specified  as 

A(1,J)  =  C(t)  ^Ao(J)  J=1,N  (4.43) 

where  A0(J)  is  defined  in  the  input  subroutine  according  to 

A0(J)  =  e»'*°,in'oJAY  j4_44j 

where  k0  is  the  wave  number  along  the  first  grid  row  and  0O  is  the  angle  the  wave  ray 
makes  with  the  x  axis.  This  offshore  condition  necessitates  that  the  water  depth  be  constant 
along  the  offshore  grid  row.  This  is  another  reason  for  choosing  the  offshore  condition  in 
the  circulation  model  to  be  a  constant  zero.  With  the  presence  of  wave  obstructions  and 
current  deflectors  such  as  jetties  it  is  possible  that  with  other  offshore  conditions  the  mean 
water  surface  elevation  could  vary  along  the  offshore  grid  row.  The  use  of  one  condition  in 
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Figure  4.2:  Start-up  function  according  to  Eq.  4.43  with  Ci=30  and  C2=3. 

the  wave  model  which  is  contrary  to  the  condition  in  the  circulation  model  can  introduce 
instabilities  in  the  overall  model.  The  interaction  of  the  two  models  is  more  sensitive  than 
either  of  the  models  individually. 

The  C(t)  term  in  equation  4.43  is  a  start-up  function  which  smoothly  goes  from  zero 
to  unity  so  as  to  bring  the  incident  wave  height  from  zero  to  H0  and  minimize  any  shock 
loading  of  the  model.  Mathematically  it  is  given  by 

C(t)  =  .5  [l  +  tanh(^-  -  C,)]  (4.45) 

where  C"i  and  C2  are  arbitrary  constants  to  adjust  the  slope  and  offset  of  the  curve.  Fig- 
ure 4.2  shows  equation  4.45  for  C*i  and  C2  equal  to  30  and  3  respectively. 

The  lateral  boundary  conditions  are  either  open  or  reflective.  The  reflective  condition 
is  expressed  as 

TJ  =  °  (4-46) 

while  the  open  condition  is  expressed  as 


3A       . 

—  =  imA 


(4.47) 


where  m  is  the  longshore  component  of  the  wavenumber  which  in  the  absence  of  diffraction 
effects  remains  constant 


m  =  k0  sin  60  =  k  sin  9 


(4.48) 
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Figure  4.3:  Waves  passing  through  a  lateral  boundary  unaffected  by  the  boundary. 

according  to  Snell's  law.  Equation  4.46  or  4.47  are  used  to  eliminate  the  amplitude  outside 
of  the  domain  from  the  tridiagonal  matrix  equation.  To  illustrate  equation  4.47  is  finite 
differenced  at  the  lower  lateral  boundary  as 


Solving  for  A(I,0) 


A(I,1)  -  A(I,0)        .     A(I,1)+A(I,0) 
AY  tm  2 


1  _  ;miY 

Mm  =       ia(i,i) 


1  +  im 


AY 


(4.49) 


(4.50) 


Thus  A(I,0)  is  eliminated  in  terms  of  A(I,1).  similarly  at  the  upper  boundary  A(I,N+1)  i 
eliminated  in  terms  of  A(I,N).  The  same  procedure  is  used  with  equation  4.46  which  yield; 
the  simple  relationships 

A(I,0)  =  A(I,1) 


and 


A(I,N+1)  =  A(I,N) 


(4.51) 


(4.52) 
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It  is  important  that  the  boundary  conditions  4.46  and  4.47  be  used  as  prescribed 
above.  The  parabolic  equation  gives  N  equations  containing  the  N+2  unknowns  A(I+1,J) 
for  J=0,N+1.  The  two  boundary  conditions  are  thus  used  to  reduce  the  number  of  un- 
knowns to  a  number  equal  to  the  number  of  equations.  If  equation  4.47  is  written  in  terms 
of  A(I+1,1)  and  A(I+1,2)  at  the  lower  boundary  and  in  terms  of  A(I+1,N-1)  and  A(I+1,N) 
at  the  upper  boundary  and  the  governing  wave  equation  is  written  for  J=2,N-1,  then  errors 
may  be  introduced  at  the  boundaries.  This  is  especially  true  for  the  case  of  a  reflective 
boundary.  Figure  4.3  shows  waves  freely  passing  through  the  lateral  boundaries  without 
any  noticeable  reflection  for  the  case  of  a  plane  wave  with  a  period  of  10  seconds  propagating 
over  constant  depth  at  an  angle  of  10  degrees  to  the  x  axis. 


CHAPTER  5 
RESULTS 


The  numerical  model  combining  the  iterative  interaction  of  the  circulation  model  and 
the  wave  model  as  described  in  the  previous  chapters  was  tested  for  several  cases  in  which 
either  analytical  or  experimental  results  are  available.  The  model  was  also  run  for  cases 
where  previous  numerical  results  are  available. 

The  two  cases  involving  analytical  solutions  are  for  set-up  induced  by  normally  incident 
waves  on  a  plane  beach  and  longshore  currents  generated  by  incident  plane  waves  at  an  angle 
on  a  plane  beach.  The  first  section  of  this  chapter  covers  set-up,  comparing  the  results 
with  the  analytical  solutions  of  Longuet-Higgins  and  Stewart  (1964)  and  the  experimental 
results  of  Bowen  et  al.  (1968)  and  Hansen  and  Svendsen  (1979).  The  second  section  covers 
longshore  currents  and  compares  results  with  the  analytical  solution  of  Longuet-Higgins 
(1970a,b)  and  experimental  results  of  Visser  (1982). 

The  model  was  run  for  the  case  of  a  rip-channel  on  a  plane  beach.  Comparisons  with  the 
numerical  results  of  Birkemeier  and  Dalrymple  (1976)  and  Ebersole  and  Dalrymple  (1979) 
are  given  in  the  third  section.  A  rip  current  was  also  generated  for  the  case  of  normal  waves 
breaking  on  a  shore  parallel  bar  with  a  gap.  No  comparisons  are  available  for  this  case. 

The  model  was  run  for  the  cases  of  shore  perpendicular  and  shore  parallel  breakwaters 
represented  in  the  model  as  infinitesimally  thin  impermeable  barriers.  For  the  shore  per- 
pendicular barrier,  the  results  are  compared  with  the  experimental  results  obtained  at  the 
Waterways  Experiment  Station  as  reported  by  Hales  (1980).  For  the  shore  parallel  barrier 
the  results  are  compared  with  the  numerical  results  of  Liu  and  Mei  (1975)  who  did  not 
include  a  wave  current  feedback  mechanism. 

In  the  final  section  the  model  was  used  to  simulate  the  laboratory  experiments  of 

63 


64 

Gourlay  (1974)  in  which  the  shoreline  in  the  lee  of  a  breakwater  was  made  to  curve  and 
meet  the  breakwater  as  occurs  with  tombolo  formation.  The  results  of  the  model  were  in 
close  agreement  with  Gourlay 's  experimental  values. 

5.1     Wave  Set-up 
Longuet-Higgins  and  Stewart  (1964)  gave  a  simple  but  elegant  derivation  for  the  wave 
induced  set-up  on  a  plane  beach  assuming  linear  shallow  water  theory.  For  the  steady  state 
case  of  normally  incident  waves  on  a  plane  beach  without  any  y  variations  and  without  any 
currents  the  depth  integrated  equations  of  motions  reduce  to 


df)   ,     1   dSxt 


r^  +  ^~5T  =  0  (51) 

In  shallow  water  for  normally  incident  waves  the  radiation  stress  term  from  linear  wave 
theory,  Eq.  (2.28),  becomes 

S„  =  E(n  +  1)  =  Z-E  =  ^P9H>  (5.2) 

A  spilling  breaker  assumption  is  made  so  that  the  wave  height  H  is  a  constant  ratio  k  to 
the  water  depth  D  =  h0  +  fj, 

H  =  K{h0  +  fj)  (53) 

Substituting  equations  (5.2)  and  (5.3)  into  Eq.  (5.1)  gives 

gg  _  3         dK?{h0  +  fjY  32/dhQ      dfj\ 

dx       ie(h0  +  fj)       Tx        -~8K  [~dx~  +  di)  <5-4) 

Solving  for  ff  gives  the  result 

gg__      f«2     dh0 

dx  l  + 1«2  dx  (5.5) 

that  the  slope  of  the  water  surface  inside  the  breaker  line  is  a  constant  times  the  slope  of 
the  beach. 

Previous  models  are  able  to  produce  the  results  of  the  analytical  expression  but  to 
date  none,  including  the  present,  have  matched  the  experimental  results.  In  figure  5.1  it 
is  seen  that  the  slope  of  the  mean  water  surface  produced  by  the  wave  set-up  is  nearly 
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Figure  5.1:  Comparison  of  the  numerical  solution  of  Vemulakonda  et  al.  (1985)  for  set-up 
with  experimental  data  of  Bowen  (1968).  Source:  Vemulakonda  et  al.  (1985) 

a  constant.   For  the  numerical  results  the  set-up  starts  at  the  breaking  point,  whereas  in 
the  experimental  results  the  water  surface  elevation  is  nearly  constant  for  some  distance 
shoreward  of  the  breaking  point  and  then  there  is  a  change  in  the  slope  of  the  water  surface 
elevation.  The  numerical  results  follow  the  analytical  solution.  This  is  an  area  that  needs 
further  investigation.  The  reason  for  the  difference  between  the  analytical  and  experimental 
results  is  that  the  analytical  formulation  is  predicated  upon  the  radiation  stress  term  being 
a  function  of  the  wave  height,  so  that  as  the  wave  height  is  reduced  after  the  breaking  point 
so  is  the  radiation  stress  term.  As  pointed  out  by  Svendsen  (1984)  the  radiation  stress  term 
remains  nearly  constant  after  breaking  for  a  distance  roughly  equal  to  five  or  six  times  the 
breaking  depth.  This  distance  is  termed  by  Svendsen  as  the  outer  or  transition  region  in 
which  there  is  a  rapid  change  of  wave  shape.  In  the  inner  region  the  wave  shape  according 
to  Svendsen  remains  nearly  constant  and  is  characterized  by  a  surface  roller  or  bore.  By 
incorporating  the  surface  roller  Svendsen  is  able  to  obtain  close  agreement  for  wave  height 
and  set-up  to  the  experimental  results  of  Hansen  and  Svendsen  (1979).  However  this  close 


agreement  is  just  for  the  inner  region.  For  the  outer  region  there  is  no  viable  solution  other 
than  that  the  radiation  stress  term  and  the  water  surface  elevation  remain  nearly  constant. 
Further  there  is  no  analytical  separation  of  the  outer  and  inner  regions  other  than  the  a 
posteri  definition  of  where  there  is  a  sharp  change  in  the  slope  of  the  water  surface  elevation. 
This  definition  is  not  useful  in  a  model  that  seeks  to  determine  the  water  surface  elevations. 
Svendsen's  work  shows  promise  and  should  be  pursued;  it  is  not  however  at  this  point 
viable  for  inclusion  in  a  numerical  circulation  model.  For  one  it  needs  to  be  expanded  for  the 
case  of  waves  at  a  angle  (this  need  will  be  discussed  further  in  the  next  section).  Another 
limitation  is  that  Svendsen's  model  is  for  monotonically  decreasing  depths,  which  is  a  severe 
restriction  (for  example  excluding  barred  profiles). 

Svendsen  conceptually  addresses  the  question  concerning  the  spacial  difference  between 
the  initiation  of  the  set-up  in  the  analytical  and  experimental  results.  The  analytical  model 
assumes  that  once  breaking  commences  energy  dissipation  also  begins.  Svendsen  states  that 
the  organized  potential  energy  of  the  wave  prior  to  breaking  is  transformed  into  forward 
momentum  in  the  outer  region.  In  keeping  with  the  definition  of  the  radiation  stresses  as 
the  flux  of  excess  momentum,  i.e.  momentum  that  would  not  be  present  in  the  absence  of 
the  wave,  it  is  easy  to  conceptualize  that  the  radiation  stress  term  should  remain  nearly 
constant  just  after  breaking.  The  problem  is  in  assuming  that  the  wave  form  in  the  breaking 
zone  is  characterized  the  same  as  outside  the  breaking  zone.  The  author,  at  present,  offers 
no  solution  other  than  to  recognize  the  existence  of  the  problem. 

In  the  present  study  results  for  the  set-up  on  a  plane  beach  with  normally  incident 
waves  differ  somewhat  slightly  from  the  analytical  solution  as  given  by  Longuet-Higgins  and 
Stewart.  This  is  because  the  present  model  does  not  use  the  spilling  breaker  assumption. 
Instead  wave  energy  is  dissipated  according  to  a  model  suggested  by  Dally  (1980)  which  is 
explained  in  section  3.3.  Thus  the  slope  of  the  mean  water  surface  is  not  a  constant  times 
the  bottom  slope.  A  result  from  the  experiments  of  Hansen  and  Svendsen  (1979)  show  this 
clearly  in  figure  5.2.  The  figure  shows  a  slight  curvature  in  the  wave  height  after  breaking 
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Figure  5.2:  Comparison  of  numerical  results  from  the  model  with  experimental  results  of 
Hansen  and  Svendsen  (1979).  Wave  set-up  and  wave  heights  for  a  2  second  wave  on  a  beach 
slope  of  1:34  and  incident  wave  height  of  7  centimeters. 

and  in  the  slope  of  the  set-up.  Although  the  problem  of  where  the  set-up  begins  is  still 
present,  the  curvature  of  the  mean  water  surface  is  the  same  for  both  cases.  It  is  also  quite 
clear  that  the  diminution  of  the  wave  height  is  not  constant  as  is  assumed  by  the  spilling 
breaker  assumption. 

In  figure  5.2  the  breaking  wave  height  in  the  experiment  is  much  greater  than  in  the 
model.  This  is  because  linear  wave  theory  cannot  adequately  predict  wave  shoaling.  Ven- 
gayil  and  Kirby  (1986)  were  able  to  duplicate  the  shoaling  recorded  in  the  experiments  of 
Hansen  and  Svendsen  (1979)  by  numerically  solving  coupled  equations  for  the  many  har- 
monic components.  Their  results  showed  that  in  the  shoaling  process  prior  to  wave  breaking, 
wave  energy  was  transmitted  to  the  higher  frequency  components.  Their  work  however  is 
not  adaptable  for  use  in  the  current  model  since  it  is  not  known  how  to  dissipate  energy 


Figure  5.3:  Non-dimensional  longshore  velocities  vresus  the  non-dimensional  distance  off- 
shore from  the  analytical  solution  of  Longuet-Higgins.  Reprinted  from  Longuet-Higgins 
(1970b),  Journal  of  Geophysical  Research,  vol.  75(33)  page  6793,  copyright  by  the  Ameri- 
can Geophysical  Union. 

from  the  various  frequency  components  once  breaking  is  initiated. 

5.2     Longshore  Currents 

Longuet-Higgins  (1970b)  gives  an  analytical  solution  for  the  longshore  current  generated 
by  obliquely  incident  waves  on  a  plane  beach.  The  solution  involves  a  spilling  breaker 
assumption  and  assumptions  concerning  the  form  of  the  bottom  friction  and  the  lateral 
mixing.  For  steady  state  conditions  the  y  equation  of  motion  reduces  to  a  balance  between 
the  gradient  of  the  Sxy  term  and  the  bottom  friction  with  the  lateral  mixing  serving  to 
distribute  the  current.  Details  of  the  solution  can  be  found  in  the  cited  paper. 

Figure  5.3  gives  the  results  of  the  analytical  solution  of  Longuet-Higgins.  In  the  figure 
the  ordinate  V  is  the  non-dimensional  velocity  v/v0,  where  v  is  the  actual  velocity  and  v0 
is  the  velocity  at  the  breaking  point  obtained  for  the  case  of  no  lateral  mixing. 

(5.6) 


5jt  i  sin  $ 

o  c 


The  abscissa  X  is  the  non-dimensional  distance  x/xh,  where  xb  is  the  distance  from  the 
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Figure  5.4:  Comparison  of  experimental  results  of  Visser  (1982)  and  numerical  result  for 
longshore  currents  due  to  a  1.1  second  wave  on  a  slope  of  1:20  for  the  deep  water  conditions 
H0  =  .065m  and  60  =  20  degrees. 

shore  to  the  breaker  line. 

Previous  models  have  reproduced  the  results  of  the  Longuet-  Higgins  solution  by  incor- 
porating all  of  the  same  assumptions  (i.e.  spilling  breakers  and  bottom  friction).  Since  the 
present  model  does  not  contain  the  same  assumptions  comparison  will  not  be  made  with  the 
analytical  results  of  Longuet-Higgins.  Instead  comparison  is  made  with  the  experimental 
results  of  Visser  (1982).  For  both  the  analytical  and  the  numerical  results  it  should  be  noted 
that  the  same  inadequacy  is  present  as  was  discussed  in  the  previous  section.  In  both  cases 
it  is  implied  that  wave  energy  dissipation  commences  as  soon  as  the  wave  begins  breaking. 
Thus  both  analytical  and  numerical  solutions  initiate  a  longshore  current  starting  at  the 
breaker  line.  It  is  reasonable  to  expect  that  the  same  holds  for  the  y  component  of  the 
onshore  flux  of  excess  momentum  as  for  the  x  component.  That  is,  that  the  Sty  term  also 
remains  nearly  constant  over  a  certain  unspecified  distance  as  the  organized  momentum  flux 
of  the  non-  breaking  wave  becomes  turbulent  momentum  upon  breaking.  This  is  illustrated 
in  figure  5.4  which  shows  a  comparison  of  the  experimental  results  of  Visser  (1982)  and  the 
model  results  for  the  same  test  conditions.  The  maximum  current  is  nearly  the  same  in 
both  the  model  and  the  experiment  but  the  maximum  is  further  offshore  in  the  model.  This 
can  be  attributed  to  the  initiation  of  a  longshore  current  in  the  model  at  the  breaker  line. 
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5.3     Waves  On  a  Rip-Channel 

Both  Birkemeier  and  Dalrymple  (1976)  and  Ebersole  and  Dalrymple  (1979)  examined 
the  case  of  waves  incident  upon  a  shore  with  the  presence  of  a  rip-channel.  The  bottom 
topography  was  given  according  to  a  formula  first  used  by  Noda  et  al.  (1974) 

Mz,y)  =  mz  jl  +  Ae-3(55)l/Ssmlo^(y-itan0)j  (5.7) 

where  £  is  the  angle  the  rip-channel  makes  with  the  x  axis  which  is,  in  their  coordinate 
system,  directed  offshore  from  the  still  water  line,  m  is  the  beach  slope,  A  is  an  amplitude 
of  the  bottom  variation  and  A  is  the  length  of  the  repeating  beach.  The  origin  of  the 
coordinates  is  at  the  intersection  of  the  still  water  line  and  the  trough  as  is  shown  in 
figure  5.5  which  gives  the  depth  contours  for  the  case  where  the  slope  is  .025  and  0  is  equal 
to  30  degrees. 

Due  to  the  restriction  of  the  present  model  that  there  be  no  y  gradients  at  the  lateral 
boundaries  the  above  formulation  is  modified  to  be 

MX)y)  =  (  mx{l  +  Ae-s(^)1/3sin10f(y-xtan^)}     |f(„  -  *tan£)|  <  * 

I  mx  l^y-rctan/?)^*  (  '  ' 

Also  shown  in  figure  5.5  are  the  contours  obtained  through  use  of  the  modified  Eq.  5.8 
for  /3  equal  to  30  degrees  and  a  .025  slope.  The  two  formulas  differ  only  in  that  the  repeating 
rip  channels  in  the  Noda  formulation  are  removed.  For  application  in  the  present  model  the 
presence  of  the  trough  in  the  farthest  offshore  grid  row  is  removed  so  that  the  depth  in  the 
offshore  grid  row  remains  constant  in  compliance  with  the  requirements  of  the  boundary 
conditions  as  discussed  in  sections  4.3  and  4.4. 

Figure  5.6  shows  the  results  obtained  by  Ebersole  and  Dalrymple  for  the  case  of  a  wave 
angle  of  30  degrees  a  wave  height  of  1.0  meters  and  a  period  of  4  seconds.  Figure  5.7  gives 
the  results  for  the  present  model.  Even  though  the  scales  of  the  two  figures  are  different  it 
is  evident  that  there  is  agreement  between  the  two.  The  length  scale  of  the  eddy  is  about 
30  meters  for  both  cases.  In  the  caption  to  the  figure  of  his  results,  Ebersole  indicated  that 
the  lateral  mixing  coefficient  may  have  been  large.   The  present  author  does  not  believe 
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Figure  5.5:  Top:  Depth  contours,  in  meters,  using  Eq.  5.7  with  the  coordinates  as  shown. 
Bottom:  Depth  contours  as  used  in  the  present  model.  Grid  spacing  is  5  meters. 
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Figure  5.6:  Velocity  vectors  from  the  results  of  Ebersole  and  Dalrymple  for  a  wave  angle 
of  30  degrees.  Scale  for  velocity  vectors  is  1  inch  =  2.87  m/sec.  Source:  Ebersole  and 
Dalrymple  (1979). 


73 


>      V     V    .» 


0.635m/s 
Maximum  Vector 


•25m 


Figure  5.7:  Velocity  vectors  obtained  from  the  present  model  using  the  same  input  condi- 
tions as  for  Ebersole  and  Dalrymple. 


E 


i  i  i  i  i  i  i  i  i 


74 


1    '    '    '    '    i    i    I    l    i    i    l    l    I    I    i    i    i    i    i    i    i 


Figure  5.8:  Mean  water  surface  elevation  contours  showing  that  the  water  surface  is  lower 
above  the  trough.  Values  indicated  are  in  millimeters. 

this  to  be  the  case.  In  the  present  study  the  lateral  mixing  coefficient  was  reduced  to  as 
low  a  value  as  2  meter2/second  without  substantially  changing  the  magnitude  of  the  eddy 
adjacent  to  the  rip  trough.  However  Ebersole  did  find  a  significant  difference  between  the 
case  of  no  lateral  mixing  and  with  lateral  mixing.  Results  for  the  present  model  for  the 
case  of  no  lateral  mixing  are  not  available  since  the  model  sometimes  exhibits  an  instability 
without  the  stabilizing  influence  of  the  lateral  mixing  terms.  Figure  5.8  shows  the  mean 
water  surface  elevation  contours  obtained  from  the  model  which  show  that  the  water  mean 
surface  is  lower  above  the  trough. 

The  agreement  between  the  results  of  the  Ebersole  model  and  the  present  model  for 
the  case  of  waves  incident  upon  a  rip-channel  should  not  be  surprising  since  diffraction  in 
this  case  is  small  compared  to  current  and  depth  refraction. 

Another  test  result  is  the  case  of  normally  incident  waves  breaking  on  a  bar  in  which 
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Figure  5.9:  Depth  contours  in  meters  for  the  case  of  a  bar  with  a  gap. 
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Figure  5.10:   Velocity  vectors  for  the  case  of  waves  over  a  bar  with  a  gap.    Dashed  lines 
indicate  depth  contours.  Grid  spacing  equal  to  10  meters. 
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there  is  a  gap  in  the  bar.  Figure  5.9  gives  the  contour  lines  for  this  test  case.  The  topography 
was  created  by  adding  an  inverted  parabola  on  to  a  plane  beach  with  a  slope  of  1:40.  The  gap 
was  created  by  suppressing  the  bar  with  a  symmetric  smooth  tangent  function.  Figure  5  10 
shows  the  velocity  vectors  resulting  from  a  6  second  wave  with  an  incident  wave  height  of  1 
meter.  The  contours  of  the  bottom  topography  are  represented  by  dashed  lines.  The  driving 
force  of  the  circulation  is  the  return  flow  through  the  gap  in  the  bar  of  the  water  that  is 
carried  over  the  bar  by  the  breaking  process.  These  results  may  not  be  totally  correct  given 
the  arguments  of  section  5.1.  If  there  were  no  gap  in  the  bar  the  waves  breaking  on  the  bar 
would  push  water  shoreward  until  sufficient  set-up  was  produced.  However  with  a  gap  in 
the  bar  lateral  flow  is  allowed  and  set-up  cannot  be  maintained. 

Kirby  (conversation  with  the  author)  has  remarked  that  he  has  observed  the  initiation 
of  a  very  strong  current  in  the  direction  of  the  waves  produced  by  waves  breaking  on  a  sand 
mound  in  a  small  wave  basin.  This  was  tested  in  the  model  using  normally  incident  waves 
on  a  slope  of  .02  with  a  shoal  described  by  the  following  equation  which  is  the  shoal  used 
in  the  Delft  experiments  of  Berkhoff  et  al.  (1982). 

•45  x  <  -5.82m 

h(x,y)=       -45  -.02(5.82  +  *)  (_^)2  +  (f)i  >  L0      (5  * 

I  •«- -02(5.82  +  x) -. 5^/1 -(jfcp -(§)»  +  . 3     fa)*  +  (|)»  <  L0 

Figure  5.11  shows  the  depth  contours  for  this  test  case.  Figure  5.12  shows  the  velocity 
vectors  obtained  for  the  case  of  a  barely  breaking  wave.  For  the  case  of  a  non-  breaking 
wave  the  velocity  vectors  plotted  to  the  same  scale  are  not  discernable.  Figure  5.13  plots 
the  max  velocity  obtained  in  the  shoal  region  versus  the  incident  wave  height.  At  a  wave 
height  of  about  10  centimeters  there  is  a  sharp  change  in  the  slope  of  the  curve  indicating 
increased  velocities  once  breaking  commences. 

5.4     Shore-Perpendicular  Breakwater 

The  model  is  tested  here  for  the  case  of  waves  interacting  with  a  shore-perpendicular 
breakwater.  The  breakwater  is  represented  as  an  infinitesimally  thin  impermeable  barrier 
located  in  between  the  J  =  JBRK  and  J  =  2BRK+i  grid  columns  and  extending  from 
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Figure  5.11:  Depth  contours  for  the  case  of  normally  incident  waves  upon  a  shoal.   Grid 
spacing  is  .05  meters.  Contour  interval  is  .02  meters. 
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Figure  5.12:  Velocity  vectors  for  the  case  of  waves  breaking  upon  a  shoal. 
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Figure  5.13:  Maximum  velocity  across  the  top  of  a  shoal  versus  incident  wave  height  for 
waves  over  a  shoal  on  a  plane  beach.  Incident  wave  height  less  than  9.1  centimeters  are 
non-breaking  waves.  Incident  wave  heights  greater  than  9.1  centimeters  produced  waves 
breaking  over  the  shoal. 
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Figure  5.14:  Position  of  the  groin  in  relation  to  the  grid  rows  and  columns. 

the  shore  boundary  to  the  I  =  \BRK  grid  row.  Figure  5.4  illustrates  the  position  of  the 
breakwater.  Mathematically,  this  thin  impermeable  barrier  behaves  as  an  internal  solid 
boundary  that  prevents  any  activity  on  one  side  of  the  barrier  from  being  transmitted  to 
the  other  side.  This  requires  modifications  to  both  wave  and  current  models. 

In  the  wave  model  an  internal  boundary  condition  is  imposed  so  that  at  the  breakwa- 
ter face  wave  energy  is  totally  reflected.  This  is  accomplished  by  imposing  the  reflection 
condition 

dA 

(5.10) 


dA 


on  both  sides  of  the  breakwater  from  grid  row  I  =  lBRK  to  the  shoreline,  or 

Mi>Jbrk)  =  A(I,Jbrk+1)  ioThlll>lBRK  (5.11) 

When  the  wave  equation  is  expressed  in  finite  difference  for  J  =  JBRK,  A{I,JBRK+1),  which 
is  on  the  other  side  of  the  barrier,  is  eliminated  by  virtue  Eq.  5.11.  Likewise  when  the  wave 
equation  is  applied  to  the  other  side  of  the  barrier  A(I,  JBRK)  is  also  eliminated.  Further 
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reflective  conditions  are  imposed  at  the  barrier  by  assuming  that  the  derivatives  of  any 
property  perpendicular  to  the  breakwater  has  a  zero  value  at  the  breakwater. 

In  the  circulation  model  the  internal  boundary  condition  of  impermeability  at  the  bar- 
rier is  imposed  by  setting  the  velocity  component  normal  to  the  barrier  equal  to 


zero. 


V{I,JBrk+i)  =  0  for  all  I  >  \BRK  (5.12) 

To  insure  that  current-induced  momentum  is  not  diffused  across  the  barrier,  subroutines 
were  introduced  in  the  computation  for  the  x  sweep  for  the  two  grid  columns  adjacent  to 
the  barrier. 

Two  sample  cases  were  run  for  comparison  purposes.  The  first  case  is  a  breakwater 
extending  400  meters  from  the  shoreline  on  a  beach  with  a  composite  slope  of  1  on  10  for 
the  first  one  hundred  meters  and  a  slope  of  1  on  100  beyond  that  point.  This  is  a  geometry 
selected  by  Liu  and  Mei  (1975).  A  deep  water  wave  height  of  1  meter  and  a  deep  water 
wave  angle  of  45  degrees  was  used  by  Liu  and  Mei.  Using  linear  shoaling  and  refraction 
this  gives  a  wave  height  of  .84  meters  and  a  wave  angle  of  31.15  degrees  as  the  wave  input 
conditions  at  the  offshore  grid  row  800  meters  from  shore  in  a  depth  of  16.95  meters. 

The  solution  method  of  Liu  and  Mei  was  to  solve  for  the  wave  field  in  the  absence  of 
currents  and  mean  water  level  changes  and  then  use  the  wave  field  as  a  constant  forcing  in  the 
solution  of  the  circulation.  By  ignoring  the  lateral  mixing  and  advective  acceleration  terms 
and  assuming  steady  state  conditions,  and  by  expressing  the  two  horizontal  components  of 
the  depth  integrated  flow  as  derivatives  of  a  stream  function  V,  the  two  horizontal  equations 
of  motion  were  reduced  to  one  Poisson  equation 

VV  =  /(*,»)  (5.13) 

where  the  /  term  contains  the  forcing  due  to  water  level  gradients,  friction  and  wave  forcing. 
This  equation  represents  a  boundary  value  problem  with  two  unknowns  ip  and  fj.  Upon 
specifying  boundary  values  an  initial  guess  for  the  water  surface  fj  is  made  and  then  Eq.  5.13 
is  solved  for  V-  With  this  new  solution  for  V  the  value  of  tj  is  updated  by  solving  steady 
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state  equations  of  motion.  This  process  is  repeated  until  convergence. 

Figures  5.15  and  5.16  show  the  results  for  Liu  and  Mei  and  the  present  study  respec- 
tively.   Figure  5.15  shows  the  contours  of  the  stream  function  whereas  figure  5.16  shows 
flow  lines  which  give  the  direction  of  the  flow  but  not  the  intensity  of  the  flow.  The  present 
model  shows  a  separation  eddy  at  the  tip  of  the  breakwater  which  is  absent  in  Liu's  results. 
Both  models  show  a  series  of  rip  cells  in  the  up  wave  side  of  the  jetty  which  result  from 
the  interaction  of  the  incident  and  reflected  waves.  There  are  not  as  many  rip  cells  in  the 
present  model  as  the  weaker  ones  farther  upstream  are  buried  by  the  stronger  longshore 
current.  Neither  the  Liu  model  nor  the  present  one  show  a  rip  current  at  the  downwave  side 
of  the  jetty.  In  figure  5.16  the  flow  lines  in  the  offshore  region  are  directed  in  the  onshore 
-  offshore  direction  which  is  different  from  the  stream  function  contour  lines  in  figure  5.15. 
This  is  because  of  a  difference  in  the  offshore  boundary  condition  in  the  two  models.   In 
the  Liu  model  the  offshore  condition  is  a  constant  xj>  which  means  that  the  flow  must  be 
tangential  to  the  offshore  boundary.  In  the  example  used  in  the  present  model  the  offshore 
condition  is  a  constant  water  surface  (fj  =  0)  which  forces  the  tangential  velocity  to  be  zero 
at  the  offshore  boundary.   An  attempt  to  change  the  offshore  condition  to  a  zero  normal 
flow  condition  did  not  accurately  conserve  water  within  the  computational  domain,  even 
though  convergence  was  achieved.  Figure  5.17  shows  the  flow  fines  resulting  from  using  no 
onshore  flow  at  the  offshore  boundary. 

The  second  case  is  to  simulate  the  experiments  conducted  at  the  Waterways  Experiment 
Station  and  reported  by  Hales  (1980).  A  jetty  perpendicular  to  the  shoreline  extended  15 
feet  out  from  the  still  water  line  on  a  plane  beach  with  a  slope  of  1  on  20.  The  plane  beach 
extended  beyond  the  jetty  tip  until  a  depth  of  one  foot  was  reached  and  further  offshore 
the  laboratory  basin  was  a  constant  depth  of  one  foot.  Waves  were  sent  at  an  angle  by 
use  of  a  moveable  wave  paddle  and  wave  guides  which  were  set  up  to  follow  the  wave  rays 
at  the  ends  of  the  wave  paddle.  Figure  5.18  shows  the  experimental  set-up  for  these  tests. 
These  wave  guides  in  essence  made  the  experiments  to  be  in  a  closed  basin  and  greatly 


83 


O    u> 
O    O 


*      V      If 
O     t*     O 

o    o 


*»  »  XM 
V  o  w 
O     o     o 


«»     fw      _       _ 

<-»       O       l«        o     f 
O      o      o       O     o 


eG=  45  r 


Figure  5.15:  Results  for  the  numerical  model  of  Liu  and  Mei  (1975).  Deepwater  wave  height 
equal  to  1  meter,  deepwater  wave  angle  equal  to  45  degrees.  Source:  Liu  and  Mei  (1975). 
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Figure  5.16:  Results  from  the  present  numerical  model  using  the  same  conditions  as  were 
used  by  Liu  and  Mei.  Offshore  boundary  condition  of  constant  water  surface  level  which 
precludes  tangential  flow  at  the  offshore  boundary. 
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Figure  5.17:  Results  from  the  present  numerical  model  using  the  same  conditions  as  were 
used  by  Liu  and  Mei.  Offshore  boundary  condition  of  no  on-offshore  flow.  Grid  spacing 
equal  to  20  meters.  6 
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Figure  5.18:   Experimental  set-up  for  the  tests  conducted  by  Hales  (1980)     Source-  Hales 
(1980). 
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Figure  5.19:   Experimental  measurements  of  the  mid-depth  averaged  velocities  obtain  by 
Hales.  Measurements  at  one  foot  intervals.  Source:  Hales  (1980). 

affected  some  of  the  circulation  patterns  away  from  the  jetty.  Figure  5.19  shows  the  average 
mid-depth  velocities  recorded  on  the  down  wave  side  of  the  jetty.  Note  the  strong  offshore 
current  next  to  the  boundary.  Since  the  boundaries  are  less  than  two  jetty  lengths  from  the 
jetty,  one  must  ask  whether  the  circulation  near  the  jetty  was  also  affected.  The  numerical 
simulation  of  the  Hales  experiments  was  performed  using  open  lateral  boundaries,  since 
to  fully  represent  the  experimental  conditions  in  which  the  lateral  boundaries  are  neither 
straight  (i.e.  curving  wave  rays)  nor  perpendicular  to  the  shoreline  boundary  would  require 
mapping  techniques  that  are  not  yet  included  in  the  capability  of  the  model. 

Figure  5.20  shows  a  comparison  of  the  experimental  and  numerical  results  obtained  for 
the  magnitude  of  the  rip  current  adjacent  to  the  down  wave  side  of  the  jetty  for  the  case  of 
a  1  second  wave  with  an  angle  of  incidence  of  30  degrees  and  a  one  foot  depth  wave  height 
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Figure  5.20:  Rip-current  velocities  adjacent  to  the  down  wave  side  of  the  jetty.  Comparison 
of  experimental  and  numerical  results. 

of  .139  ft.  One  would  expect  the  depth  averaged  velocity  to  be  somewhere  between  the 
mid-depth  velocity  and  the  bottom  velocity.  The  model  appears  to  slightly  underpredict 
the  rip  current.  In  figure  5.19  a  circulation  cell  in  the  lee  of  the  groin  is  barely  discernable. 
In  figure  5.21  the  circulation  cell  is  very  obvious. 

5-5  Shore  Parallel  Emerged  Breakwater 
In  order  to  solve  for  the  wave  and  current  field  resulting  from  the  interaction  of  waves 
and  an  emerged  breakwater  several  boundary  conditions  at  the  breakwater  are  incorporated 
into  the  model  in  addition  to  the  boundary  conditions  already  discussed  in  sections  4.3 
and  4.4.  The  breakwater  is  represented  as  being  infinitesimally  thin  and  impermeable. 
Figure  5.22  shows  the  position  of  the  breakwater  on  the  shoreward  side  of  the  I  =  lBRK  grid 
row.  The  impermeability  condition  implies  that  neither  waves  nor  currents  pass  through 
the  breakwater.  The  no  current  condition  is  easily  imposed  in  the  circulation  model  by 
making  U(lBRK+1,  J)  equal  to  zero  for  the  entire  length  of  the  breakwater.  In  addition  the 
circulation  model  is  modified  so  that  momentum  diffusion  and  derivatives  do  not  take  place 
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Figure  5.21:  Velocity  vectors  of  the  currents  obtained  numerically  for  the  same  conditions 
as  for  the  Hales  experiment.  Grid  spacing  equal  to  one  foot. 
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Figure  5.22:  Definition  sketch  of  grids  for  the  case  of  a  shore  parallel  breakwater, 
across  the  barrier. 

In  the  wave  model  an  adaptation  is  made  so  as  to  prescribe  the  condition  that  the 
wave  height  is  zero  on  the  lee  side  of  the  breakwater.  This  condition  on  the  lee  side  of  the 
breakwater  is  physically  not  realistic  since  common  observation  reveals  that  waves  actually 
propagate  along  the  lee  side  of  a  breakwater.  The  parabolic  model  however  marches  waves 
in  but  one  direction.    This  assumption  however  will  give  good  results  several  grids  away 
from  the  breakwater.   Figure  5.23  shows  the  crests  of  the  diffracted  waves  in  the  lee  of  a 
semi-infinite  breakwater  for  a  region  of  constant  depth.    This  figure  shows  that  the  zero 
wave  height  condition  on  the  lee  side  of  the  breakwater  does  not  adequately  represent  the 
waves  along  the  breakwater  or  allow  wave  energy  to  properly  propagate  in  the  lee  of  the 
breakwater;  however,  important  information  can  be  obtained  through  use  of  the  marching 
parabolic  wave  equation  method.    To  set  the  wave  height  to  zero  on  the  lee  side  of  the 
breakwater  the  ACOEF  subroutine  is  modified  so  as  to  advance  the  wave  one  half  delta  x 
from  the  center  of  the  I  =  BRK  grid  row  to  the  position  of  the  breakwater  at  the  edge  of 
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Figure  5.23:  Wave  crests  in  the  lee  of  a  semi-infinite  breakwater  in  a  region  of  constant 
depth. 

the  grid  row.  The  amplitude  is  then  set  to  zero  at  the  breakwater  and  remains  the  same 
elsewhere.  The  wave  is  then  advanced  another  one  half  delta  x  so  as  to  be  back  at  the  grid 
centers  and  then  is  marched  grid  row  by  grid  row  to  the  shoreline.  With  the  exception  of 
the  two  one  half  delta  x  steps  the  wave  model  is  otherwise  unchanged. 

A  test  case  was  arbitrarily  chosen  of  a  breakwater  100  meters  long  100  meters  from 
the  shoreline  on  a  plane  beach  with  a  1  on  40  slope.  The  wave  conditions  again  were 
arbitrarily  chosen  to  be  a  sue  second  period  and  a  wave  height  of  1  meter  on  the  offshore 
grid  row.  Figure  5.24  shows  the  velocity  vectors  obtained  for  the  case  of  normal  incidence 
and  for  the  cases  of  increasing  angle  of  incidence.  For  normal  incidence  two  circulation  cells 
and  a  well  defined  rip  current  in  the  lee  of  the  breakwater  result.  The  rip  current  is  fed 
by  the  longshore  currents  which  result  from  the  difference  in  the  set-up  in  the  lee  of  the 
breakwater  and  the  region  exposed  to  direct  waves.  As  the  angle  of  incidence  increases  a 
longshore  current  develops  and  the  circulation  pattern  becomes  skewed.  The  rip  current 
is  increasingly  skewed  in  the  direction  of  the  longshore  current  with  increasing  strength  of 
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the  longshore  current  and  the  upstream  circulation  cell  becomes  increasingly  weaker.  The 
other  cell  remains  but  is  shifted  down  current. 

The  emerged  breakwater  program  was  also  run  for  the  case  of  a  shore  parallel  breakwater 
700  meters  in  length,  350  meters  offshore  on  a  1  on  50  slope.  This  is  a  case  tested  by  Liu 
and  Mei  and  represents  the  condition  for  the  Santa  Monica  offshore  breakwater.  The  wave 
conditions  were  a  10  second  wave  of  normal  incidence  with  a  deep  water  wave  height  of  1 
meter.  Figure  5.25  shows  the  symmetrical  stream  function  contour  lines  obtained  from  the 
Liu  model  for  half  of  the  region  in  the  lee  of  the  breakwater.  Figure  5.26  shows  the  flow  lines 
obtained  from  the  present  model  and  figure  5.27  shows  the  velocity  vectors  for  this  case. 
Figure  5.28  shows  the  contour  fines  of  the  mean  water  surface  or  set-up.  There  is  obviously 
greater  set-up  in  the  exposed  region  than  in  the  lee  of  the  structure  but  no  well  defined  rip 
current  is  seen  in  figure  5.27  or  in  the  results  of  Liu  shown  in  figure  5.25.  Both  models  show 
a  strong  longshore  current  at  the  shoreline  directly  behind  the  end  of  the  breakwater,  but 
this  current  does  not  extend  to  the  center  of  the  breakwater  since  there  is  a  distance  over 
300  meters  in  which  the  intense  longshore  current  can  be  diffused  and  directed  offshore. 
5.6     Comparison  with  Experimental  Results  of  Gourlay 

Gourlay  (1974)  reported  results  of  laboratory  experiments  in  which  currents  were  gen- 
erated by  the  diffracted  waves  in  the  lee  of  a  breakwater.  His  experiment  was  designed  to 
demonstrate  that  longshore  variation  of  the  wave  height  could  generate  longshore  currents. 
The  experimental  wave  basin  was  thus  constructed  so  that  the  shoreline  was  everywhere  par- 
allel with  the  diffracted  wave  crests.  Figure  5.29  shows  the  laboratory  set-  up  for  Gourlay 's 
experiments.  A  1  on  10  concrete  beach  is  parallel  to  the  incoming  wave  crests  in  the  ex- 
posed zone  and  in  the  shadow  zone  the  slope  is  also  1  to  10  but  the  beach  is  curved  with 
a  constant  radius  centered  on  the  breakwater  tip.    Gourlay's  reasoning  in  choosing  this 
arrangement  was  that  "since  the  diffracted  wave  crests  have  theoretical  curves  which  can 
be  approximated  by  a  circle  the  influence  of  waves  breaking  at  an  angle  to  the  beach  was 
almost  completely  eliminated."  (page  1978) 
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Figure  5.25:  Stream  function  contour  lines  obtained  from  the  Liu  model  for  a  10  second  1 
meter  wave  normally  approaching  a  700  meter  long  breakwater  350  meters  offshore  on  a  1 
on  50  slope.  Source:  Liu  and  Mei  (1975). 
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Figure  5.26:    Flow  lines  obtained  from  the  present  model  for  a  10  second  1  meter  wave 
normally  approaching  a  700  meter  long  breakwater  350  meters  offshore  on  a  1  on  50  slope. 
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Figure  5.27:  Velocity  vectors  obtained  from  the  present  model  for  a  10  second  1  meter  wave 
normally  approaching  a  700  meter  long  breakwater  350  meters  offshore  on  a  1  on  50  slope. 
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Figure  5.28:  Set-up  contour  lines  obtained  from  the  present  model  for  a  10  second  1  meter 
wave  normally  approaching  a  700  meter  long  breakwater  350  meters  offshore  on  a  1  on  50 
slope. 

An  experiment  of  Gourlay  was  numerically  duplicated  by  making  several  adaptations 
to  the  overall  scheme.  Figure  5.30  shows  the  numerical  computational  grids  used  in  the  pro- 
gram. Rather  than  model  the  entire  wave  basin  it  was  decided  for  computational  economy 
to  assume  the  wave  paddle  location  to  be  situated  3  meters  offshore  of  the  breakwater  and 
that  the  wave  paddle  produced  waves  only  in  half  of  the  domain  seaward  of  the  breakwater. 
The  offshore  boundary  condition  in  the  circulation  model  was  a  no-flow  condition  through 
the  wave  paddle.  This  condition  may  slightly  influence  the  results  obtained  for  the  mean 
water  levels  since  the  computer  model  will  not  have  as  great  an  off  shore  reservoir  from 
which  the  wave  set  up  must  come  as  in  the  physical  experiment.    An  adaptation  to  the 
flooding  routine  was  made  to  allow  for  variable  flooding.    This  meant  that  the  last  grid 
of  each  column  could  vary  independently  of  the  neighboring  grid  columns.    This  change 
necessitated  the  incorporation  of  many  decision  statements  within  the  circulation  model  so 
that  there  would  be  no  flow  between  wet  grids  and  dry  grids  unless  the  water  elevation  was 
high  enough  to  initiate  flooding.  At  present  flooding  is  allowed  only  in  the  x  or  on-offshore 
direction. 
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Figure  5.29:  Physical  layout  for  the  experiment  of  Gourlay.  Source:  Gourlay  (1974),  Pro- 
ceedings 14  Coastal  Engineering  Conference,  copyright  by  the  American  Society  of  Civil 
Engineers,  reprinted  with  permission. 


98 


X 


— 


— 


Figure  5.30:  Computational  domain  for  the  experiment  of  Gourlay.  Grid 
.1  meter. 
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In  the  wave  model  part  of  the  program,  the  problem  of  wet  and  dry  grids  is  handled 
using  a  technique  proposed  by  Kirby  and  Dalrymple  (1986)  by  treating  the  dry  grids  as 
though  they  have  a  very  small  (say  one  millimeter)  depth  of  water.  This  technique  allows  for 
keeping  a  uniform  computational  domain  for  the  wave  model  while  avoiding  any  numerical 
problems  such  as  dividing  by  zero  depths. 

Gourlay's  experiment  used  a  wave  height  of  9.1  centimeters  with  a  wave  period  of  1.5 
seconds.  The  depth  in  the  offshore  constant  depth  region  was  20  centimeters.  For  the 
numerical  computation  a  grid  size  of  .1  meters  was  used.  This  created  a  60  by  63  grid 
starting  domain  which  grew  to  a  64  by  63  grid  domain  once  set  up  was  established. 

The  numerical  results  show  close  agreement  with  the  experiment  results.  Figure  5.31 
shows  the  experimental  results  for  the  mean  water  surface  contours  resulting  from  the  wave 
set-up  and  the  corresponding  results  obtained  from  the  numerical  model.  The  agreement 
appears  to  be  quite  good.  There  is  however  an  area  of  slight  differences  where  the  shoreline 
connects  to  the  breakwater.  The  computer  model  fails  to  show  any  set-up  in  contrast  to  the 
experimental  results  which  show  a  slight  set-up  in  this  region.  An  explanation  for  this  is 
that  the  parabolic  approximation  used  in  the  wave  model  does  not  propagate  waves  along 
the  breakwater  as  occurs  in  the  physical  experiment.  This  in  turn  may  also  help  explain 
the  absence  of  the  small  stagnation  eddy  shown  in  figure  5.32. 

Figure  5.32  shows  the  experimental  results  for  the  circulation.  The  figure  shows  both 
the  contours  for  the  magnitude  of  the  velocity  and  some  streamlines  indicating  the  direction 
of  the  flow.   Figure  5.33  shows  the  stream  lines  obtained  from  the  numerical  values  while 
figure  5.34  shows  the  contour  lines  for  the  magnitude  of  the  velocities  obtained  numerically. 
Again  good  overall  agreement  is  in  evidence.  The  rip  current  in  the  exposed  area  in  par- 
ticular is  very  well  shown  in  the  numerical  results.  The  differences  are  that,  as  mentioned 
earlier,  there  is  no  stagnation  eddy  shown  in  figure  5.33  and  the  center  of  the  primary  eddy 
is  slightly  farther  behind  the  breakwater  in  the  physical  experiment  than  in  the  model  re- 
sults. Similarly  as  seen  in  figure  5.34  the  maximum  currents  obtained  in  the  model  are  not 
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Figure  5.31:    Comparison  of  experimental  and  numerical  results  for  the  experiment  of 
Gourlay:    Wave  Set-up  Contours.     Top  figure  is  from  Gourlay  (1974)  Proceedings  14th 
Coastal  Engineering  Conference,  copyright  by  the  American  Society  of  Civil  Engineers 
reprinted  with  permission.  The  bottom  figure  is  the  computer  results. 
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Figure  5.32:  Results  from  the  experiment  of  Gourlay:  Contours  of  current  velocities  and 
streamlines.    Source:   Gourlay  (1974),  Proceedings  14th  Coastal  Engineering  Conference 
copyright  by  the  American  Society  of  Civil  Engineers,  reprinted  with  permission 
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Figure  5.33:  Results  from  the  numerical  model:  Streamlines  of  the  flow. 
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Figure  5.34:  Results  from  the  numerical  model:  Contour  lines  of  the  current  velocities. 

as  far  behind  the  breakwater  as  they  were  in  the  experiment.  Both  of  these  spacial  dispar- 
ities between  numerical  and  physical  results  may  result  from  the  failure  of  the  parabolic 
approximation  to  adequately  propagate  wave  energy  fully  into  the  lee  of  the  breakwater. 
In  comparing  the  velocities  obtained  in  the  model  to  the  experimentally  measured  values  it 
should  be  noted  that  Gourlay's  results  are  for  surface  velocities  whereas  the  model  results 
are  for  depth-averaged  currents.  This  may  be  a  factor  in  explaining  that  the  maximum  in 
figure  5.32  for  the  experimental  results  is  larger  than  the  maximum  in  figure  5.34  for  the 
numerical  results. 

The  differences  in  the  experimental  and  numerical  results  which  were  mentioned  above 
should  be  kept  in  perspective.  Overall  the  agreement  is  good  both  qualitatively  and  quan- 
titatively. 


CHAPTER  6 
CONCLUSIONS  AND  RECOMMENDATIONS  FOR  FURTHER  STUDY 


A  numerical  solution  for  the  circulation  and  waves  in  coastal  regions  has  been  presented. 
The  solution  method  uses  an  alternating  direction  implicit  method  of  solving  the  depth 
integrated  equations  of  momentum  and  continuity.  Wave  forcing  terms  are  obtained  through 
a  parabolic  approximation  to  the  linear  mild  slope  equation  for  wave  propagation  over 
varying  depths  and  currents.  The  numerical  solution  involves  an  iteration  between  the 
parabolic  wave  solution  and  the  ADI  circulation  solution. 

The  model  was  run  for  several  test  cases.  These  include  cases  with  varying  topogra- 
phy, such  as  rip  channels  and  bar  gaps.  The  model  was  also  tested  with  shore  parallel  and 
shore  perpendicular  breakwaters  represented  as  infinitesimally  thin  impermeable  barriers. 
Comparisons  were  made  with  analytical  and  experimental  results.  These  comparisons  were 
favorable.  In  particular  comparison  of  the  results  from  the  model  in  simulating  the  experi- 
ments of  Gourlay  (1974)  show  the  predictive  capability  of  the  model.  The  model  was  also 
run  for  an  idealized  shore  parallel  breakwater  varying  the  aspect  ratio  of  the  length  of  the 
breakwater  to  the  distance  offshore  of  the  breakwater  to  determine  the  relative  importance 
of  this  parameter. 

The  primary  conclusion  of  this  study  is  that  it  is  possible  to  get  good  predictive  re- 
sults using  a  numerical  model  that  incorporates  some  simplifications  of  the  physical  reality. 
In  the  present  study  the  three-dimensional  flows  that  are  encountered  in  nature  are  rep- 
resented by  two-dimensional  depth  averaged  velocities.  Simplifying  assumptions  are  also 
made  concerning  bottom  friction,  momentum  transfer  or  mixing,  wave  breaking  and  wave 
energy  dissipation.  For  each  of  these  topics  it  is  possible  to  argue  that  nature  is  greatly  dif- 
ferent and  more  complex  than  the  assumptions  embodied  in  this  report.  Yet  the  results  (in 
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particular  the  comparison  to  the  laboratory  experiments  of  Gourlay)  confirm  the  predictive 
capability  of  the  numerical  model.  The  simplifications  mentioned  above  make  this  possi- 
ble, for  there  must  be  a  mathematical  statement  of  the  physical  process  before  numerical 
modelling  of  the  process  can  even  be  attempted. 

In  making  recommendations  for  further  study  we  should  distinguish  between  short  range 
objectives  and  the  possibilities  that  may  become  open  to  numerical  modelers  by  advances 
in  the  computer  hardware  field.  It  is  not  difficult  to  envision  within  a  few  decades  or  less, 
super  computers  with  thousands  of  parallel  processors.  Such  advances  will  make  many  of  the 
present  computational  methods  obsolete.  For  example  the  parabolic  approximation  is  made 
solely  to  put  the  wave  equation  in  a  solvable  form  with  respect  to  the  capabilities  of  present 
day  computers.  With  sufficient  parallel  processors  it  is  conceivable  that  it  may  be  practical 
to  solve  the  full  elliptic  equation  explicitly.  We  shall  not  dwell  further  on  speculation  about 
future  generations  of  computers.  In  the  near  term,  the  numerical  model  presented  in  this 
report  can  be  refined  and  the  applications  can  be  extended.  The  model  is  constructed  in 
modular  form  so  that  any  new  formulation  for  any  of  the  component  processes  can  easily 
be  incorporated  into  the  model.  The  model  is  thus  open  for  continual  improvement. 

6.1     Improved  Wave  Formulations 

The  model  can  be  improved  by  better  wave  modeling  and  by  a  more  correct  formulation 
of  the  radiation  stress  terms  once  breaking  has  commenced.  The  formulation  used  in  this 
model  and  previous  models  for  the  radiation  stresses  initiates  set-up  and  longshore  currents 
at  the  breaker  line  whereas  as  explained  in  section  5.1  experiments  show  otherwise. 

An  improved  wave  model  will  assign  a  non-zero  value  of  the  amplitude  in  the  lee  of 
a  shore  parallel  breakwater.  Research  should  be  done  to  establish  a  method  of  doing  this 
so  as  to  enable  more  wave  energy  to  propagate  into  the  lee  of  the  structure.  Dalrymple 
and  Kirby  (1988)  have  found  a  method  to  determine  a  non-zero  wave  height  in  the  lee  of  a 
breakwater.  However  their  solution  is  in  the  form  of  a  Fourier  integral  which  is  not  suited 
to  the  complex  amplitude  form  required  in  the  present  model. 
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6.2     Including  Wave  Reflection 

The  author  has  made  unsuccessful  attempts  to  include  wave  reflection  in  the  model. 
Kirby  (1986)  used  coupled  equations  governing  the  interaction  of  a  forward  and  a  backward 
(reflected)  propagating  wave  that  are  coupled  through  a  slope  interaction.  He  was  able  to 
obtain  the  reflection  and  transmission  coefficients  over  submerged  breakwaters  that  agree 
well  with  the  experiments  of  Seelig  (1980).  Once  the  arrays  of  the  complex  amplitudes  of 
the  two  waves  is  determined,  the  radiation  stress  terms  can  be  determined  from  equations 
that  have  been  derived  following  the  same  procedure  as  used  by  Mei  (1972)  but  defining  4> 
as  being  comprised  of  two  components.  In  fact  this  procedure  can  also  be  used  to  obtain 
the  radiation  stress  terms  for  any  number  of  components.  The  only  limit  being  the  amount 
of  computation  required  as  the  number  of  components  increases. 

Smith  (1987)  following  the  work  of  Miles  (1967)  for  wave  reflection  and  transmission 
over  a  step  derived  equations  for  the  complex  reflection  and  transmission  from  a  submerged 
plane  barrier.  An  attempt  was  made  using  this  formulation  to  obtain  the  circulation  and 
wave  pattern  from  a  submerged  breakwater  represented  as  a  plane  impermeable  barrier. 
However  problems  of  numerical  stability  were  encountered.  It  is  suspected  that  the  sharp 
gradients  between  the  region  with  the  presence  of  the  reflected  wave  and  regions  of  its 
absence  adjacent  to  the  barrier  produced  unrealistic  gradients  in  the  radiation  stress  terms. 
An  attempt  was  made  to  include  the  reflected  wave  for  the  case  of  an  emerged  plane  barrier 
with  the  same  unstable  result.  This  is  a  topic  for  further  investigation. 

63     Extending  the  Model  to  Include  Sediment  Transport 

The  model  as  it  currently  stands  can  be  extended  to  include  sediment  transport.  The 
wave  model  can  be  slightly  amended  to  retain  wave  energy  dissipation  as  an  array.  If  a 
formulation  of  sediment  entrapment  can  be  made  in  terms  of  the  wave  energy  dissipation 
and  the  bottom  shear  stress,  everything  else  is  in  place  to  accomplish  the  sediment  transport 
model.  A  simple  numerical  mechanism  would  be  to  assume  steady  state  conditions  for 
certain  time  intervals  and  solve  for  the  wave  and  current  field  and  the  sediment  concentration 
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as  a  function  of  known  variables  and  then  solve  a  separate  model  for  sediment  balance  and 
then  update  the  depths  and  proceed  iterating  between  the  sediment  balance  model  and  the 
present  wave-current  model  marching  in  time.  Since  it  has  already  been  demonstrated  that 
the  wave-  current  model  gives  adequate  results  in  the  lee  of  a  shore  parallel  breakwater,  a 
good  test  for  such  a  sediment  transport  model  wold  be  to  see  if  it  could  yield  a  tombolo 
formation  in  the  lee  of  a  shore  parallel  breakwater. 


APPENDIX  A 
DERIVATION  OF  THE  DEPTH  AVERAGED  EQUATIONS 


In  this  appendix  a  derivation  of  the  governing  equations  for  the  circulation  model  is 
presented.  This  appendix  closely  follows  the  work  of  Ebersole  and  Dalrymple  (1979). 

Before  proceeding  with  the  derivations  of  the  depth  integrated  time  averaged  equations 
of  motion  and  the  continuity  equation  the  Leibnitz  Rule  for  taking  the  derivative  of  an 
integral  is  stated  since  it  is  invoked  throughout  the  derivation. 

The  Leibnitz  Rule  can  also  be  used  in  reverse  to  take  derivatives  outside  of  an  integral. 
The  kinematic  boundary  conditions  are  also  reproduced  since  the  use  of  the  Leibnitz  Rule 
will  yield  bottom  and  surface  terms  that  will  be  eliminated  through  use  of  these  boundary 
conditions.  On  the  free  surface  the  boundary  condition  is 

dr\  dr\  drj 

and  at  the  bottom  the  boundary  condition 

dh0  dh0  dh0 

In  the  circulation  model  it  is  assumed  that  the  still  water  depth  does  not  vary  with  time  so 
the  first  term  of  Eq.  (A.3)  is  dropped  (actually  for  time  varying  depth  essentially  the  same 
equations  are  obtained). 

The  continuity  equation  (or  the  conservation  of  mass  equation  as  it  is  also  called)  in 
three  dimensions  is 

gg  |   gtgO   ,   d(pv)      djpw) 
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Integrating  over  the  depth  from  z  =  -k0  to  z  =  q  and  using  the  Leibnitz  Rule  to  take 
derivatives  outside  the  integral,  the  continuity  equation  becomes 


,     x   drj  dh0      d    f> 


fin  Ah  ^ 

-(H^-(H-A<,^  +  (^)n-(£uO_^    =    0  (A.5) 

The  terms  with  an  overbrace  are  p  times  the  free  surface  kinematic  boundary  condition 
and  the  terms  with  an  underbrace  are  p  times  the  bottom  boundary  condition.  Eliminating 
both  of  these  from  Eq.  (A.5)  gives 

Ft  U  pdz  +  rx  U  pudz+rv  U pvd2  =  0  (A6) 

Equation  (A.6)  is  now  time  averaged.    Time  averaging  means  averaging  over  a  wave 
period.  The  time  average  of  a  term  is  denoted  by  an  overbar.  Thus 

J=fIoFdt  <A-7) 

where  F  represents  any  term  and  T  is  the  wave  period.  Letting  u  and  v  be  comprised  of  a 
mean  flow,  a  wave  induced  flow,  and  an  arbitrary  or  random  turbulence  component 

u  =  u  +  ti  +  u'  (A.8) 

v  =  v  +  v  +  v'  (A.9) 

The  turbulent  fluctuations  are  assumed  to  be  of  a  high  frequency  relative  to  the  wave 
frequency  so  that  their  time  averages  over  the  wave  period  are  zero.  By  substituting  the 
above  expressions  for  u  and  v  into  Eq.  (A.6)  and  taking  time  averages  the  continuity  equation 
is 

Jl  {P(h0  +  fj)}  +  —  {pu(h0  +  fj)}  +  —  j*    ptdz+  (A.10) 


J-y  {*>(*.  +  *»  +  d-yf-hopidz  =  °  (A11) 
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The  time  averages  of  the  vertically  integrated  wave  induced  velocities  u  and  v  are  not  zero. 
These  quantities  are  known  as  the  wave  induced  mass  flux. 
Defining 

U  =  fi  +  u  and  V  =  v  +  v  (A. 12) 

where 


u  = 


(A.13) 


the  time  averaged  depth  integrated  continuity  equation  is 

Tt  +kU(K  +  f)  +  ^(*.  +  *)-0  (A.14) 

The  assumption  that  the  density  is  constant  in  both  time  and  space  has  used  in  obtaining 
both  Eq.  (A. 11)  and  Eq.  (A.14).  The  assumption  of  a  constant  bottom  has  also  been  used 
in  obtaining  (A.14).  If  the  bottom  is  not  constant  in  time  a  ^f  term  would  be  included  in 
Eq.  (A.14). 

The  x  and  y  momentum  equations  are  handled  in  the  same  manner  as  the  continuity 
equation  above.  The  x  and  y  direction  momentum  equations  are 

du      dtS       duv      duw  _  1  dp       1    drt *       3rvx       drMt , 

dt^  dx  +  dv  +  dz  ~  0di  +  ~o{~dx-  +  -a^  +  -fi7i  (A15) 


^      ^      M       dw       Idp      Idr^      dry        dr, 

dt*   dx   *  dy  +    dz         pdy  +  p{  dx    +1y~+~dT}  (A16) 

The  same  procedure  is  followed  for  both  the  x  and  y  equations.  Only  Eq.  (A.  15)  will  be 
dealt  with  here.  Integrating  the  left  hand  side  over  the  depth  and  applying  the  Leibnitz 
Rule  to  take  the  derivative  outside  the  integral  gives 

(LHS)   "   dll-Cdz~u«Tt-uj^jl+ 

J -h0  ''dx        -*•  dx 

(uw),,  -  {uw).ho  (A.17) 


a_ 

dx 
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The  terms  in  Eq.  (A.17)  that  are  designated  by  an  overbrace  are  u„  times  the  kinematic 
free  surface  condition  and  those  designated  by  an  underbrace  are  u_fco  times  the  bottom 
boundary  condition.  These  are  both  dropped  leaving 

(LHS)  =  §iJ"ho"**+£j"h/d*+§-J"hUvdz  (A.18) 

Substituting  the  definitions  given  for  u  and  v  in  Eqn.  (A.8)  and  (A.9)  and  recognizing  that 
the  turbulent  fluctuating  components  are  of  a  high  frequency  relative  to  the  wave  frequency 
so  that  the  time  averages  of  integrals  involving  the  product  of  a  turbulent  component  and 
a  wave  induced  component  are  zero  gives 

LHs "  HZ. s  * + mil.  •  * + j-XK  •** + iLLf"*    <ai9> 

+  Yxihy^dz  +  aiU  2™dz  +  aiU  ™dz  +  aiU  ™dz         (A20) 

B  ~H>    "  d  ~pi  d  ~n 

+  TyL.  U€  ^  +  TvU  » dz  +  TyU  «'*' dz  (A21) 

using  the  definitions  of  Eqs.  (A.12)  and  (A.13)  this  becomes 


LHS    =    |^  +  ,,  +  A^  +  ,-)+^^_^[_i_{/_:/i<(2} 


+hm^ « +fvi:. « *  -  £  y^ij^z? * 


(A.22) 


where  the  two  terms  designated  by  underbraces  result  from  being  added  to  complete  the 
U2  and  UV  terms. 

The  right  hand  side  after  integrating  over  the  depth  assuming  a  constant  density  and 
assuming  an  inviscid  fluid  such  that  no  horizontal  viscous  stresses  exist,  i.e.  r«  and  ryx  are 
zero,  is 

p  J-ho  dx      T  p  J.Ke   dz  aZ  ^2Z) 

Using  the  Leibnitz  Rule  and  averaging  over  a  wave  period  the  right  hand  side  becomes 


p  dx].ho  P      T  pP"dx  +  pP~h°  dx  +  pT""       pT"-»c 


(A.24) 


Ill 

r^"  and  r„_ko  are  the  time  averaged  surface  and  bottom  shear  stresses.  Assuming  that 
the  free  surface  pressure  p„  is  zero  (the  dynamic  free  surface  boundary  condition)  the  right 
hand  side  is 


RHS"  -pTxUpdz  +  -/-»•  *T+  ~^~  -^  (A.25) 

The  mean  pressure  at  the  bottom  can  be  expressed  as  the  sum  of  the  wave  induced  or 
dynamic  pressure  and  the  hydrostatic  pressure 

P^ho  =  Pdyn.ko  +  pg{h0  +  fj)  (A.26) 

so  that  the  product  p^^fe  can  be  manipulated  to  read 

v=s^k  -  \^^d~t + \ih{9{ho + n),}  ~  ^ +^Yx       (A-27> 

Terms  from  both  the  right  and  left  hand  sides  are  combined  as  the  radiation  stresses. 
Defining 


s" ~ I -J* + p"2) dz " 2pg(h° + ^ ~ pjir+^itf-  p<ldz^      (A-28) 

S,v  =  f\utdz  -  -^J^J^Z  (A.29) 

Sw  =  J_hy+f°2) dz  -  2P9^ho  +  *>?  ~  nr+  -) ij    Pidz)2       (A-3°) 

and  making  the  assumptions  that  the  product  of  the  dynamic  pressure  at  the  bottom  and 
the  bottom  slope  is  negligible,  that  the  &f2ko(u')*  dz  term  can  be  ignored  and  that  the 
lateral  friction  rt  caused  by  momentum  fluxes  due  to  turbulent  fluctuations  is  independent 
of  depth  and  that  r,  is  defined  by  r,  =  -pu'v',  the  x  direction  equation  of  momentum  is 
written  as 

*U(ha  +  If)  +  §-U\h0  +  fj)  +  ±UV(h0  +  fj)  +  l*p*  +  l*fJ3L 
at  ox  dy  p   dx        p    dy 

,  _.dfj      {h0  +  fj)dTt       1  1 

+,(*.  +  „)_  +  i-^-i^  -  -r„  +  -*.  =  0  (A.31) 

By  subtracting  the  depth  integrated  continuity  equation  and  dividing  by  the  total  depth, 
the  final  form  of  the  x  momentum  equation  is 

^  +  u^  +  vJy-  +  pJ^TT)[-aT  +  -af}  (A32) 
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+gdx  +  pd^-  7U^TT)T"  +  7{KTW)n*  =  °  (A-33) 

The  same  procedure  gives  as  the  depth  integrated  time  averaged  y  momentum  equation 

^  +  u—  +  VdV  +         ?         \dS">  ^  dS*v1 
dt+UdX+V^  +  rt^r)[-d?  +  -e?\  (A.34) 

+gdy +-p*;-  pv^T)r-« + pjh^nrf* = °  (A-35) 


APPENDIX  B 
DERIVATION  OF  THE  RADIATION  STRESS  TERMS 


The  basis  of  the  derivation  is  to  use  the  definition  of  the  velocity  potential  to  substitute 
for  the  wave  induced  velocities  in  the  definition  for  the  radiation  stresses. 

d<j> 


«i- 


dxp 


,       8<f> 

w=Tz 


where  (p  is  expressed  in  terms  of  the  general  complex  amplitude  function  B 


as 


where  ex.  represents  the  complex  conjugate  term.  Also  r}'  is  expressed  as 

rfml  (Be"*-*  +  c.c) 

where  the  prime  is  used  to  denote  wave  induced  quantities. 

For  convenience  Eq.  (2.50)  of  Section  2.2  is  reproduced  below. 

8+  =  pf^  u^dz  +  SaP  [J%dz  -P{(h0  +  fi?  +  f  W] 
An  expression  for  p  is  obtained  in  Section  2.2  and  is  reproduced  bel 


ow 


P[x,y,z)  =  pg(r)-z)  +  p        —L—dz-p(w>y 

Jm      uX( 

Substituting  Eq.  (B.5)  into  Eq.  (B.4)  gives 


Sap    =    P  u'au'0dz  +  6af) 


(B.l) 


(B.2) 


(B.3) 


(B.4) 


(B.5) 


-/_fc/K)2^-f(/»o  +  ^  +  fw 


(B.6) 
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The  terms  of  Eq.  (B.6)  are  now  evaluated  term  by  term. 
Term  1: 


a»  =  I  (_*9dB  cosh k(h0  +  z)      ■  ,      ig3B' cosh  *(*.  +  *)  M\ 
"      2\    aflia       coshjfcfc  ^  a  dxa        c^sh**       *     J  (BJ) 


i  mi  _  g2  cosh2  k(h0  +  z) 
4<r2  cosh2  kh 


«««/>  = 


dxa  dxp       dx*  dip      dxa  dxp  dxa  dxfi  * 


(B.8) 


pW    =     P9  cosh2  k{h0  +  z) 
a  fi  2Jfcsinh2A:/i 

pg  cosh2  k(h0  +  z) 
k  sinh  2kh 


dB  dB*      dB"  dB 


dxa  dip       dxa  dip 
dB  dB 


SR 


dxa  dip 


; 


(B.9) 


/      pu'au'pdz    = 

J  —h0 


ksinh2kh 
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ksmh.2kh 


'(  dB  dB*\\   ft  , 

=  nm(dBBB*\  i  /       abfc  \ 

4      \dxadx{,J  k2  V        sinh  2**/ 


(B.10) 


The  second  term: 


y_fc  M*? -*)<**=  pg(^-y) 


—  w — « — 


-Ac 


(B.ll) 


The  third  term: 


)eiu,t 


u'     B     1  /    «g  aB  cosh  fc(ft0  +  z)  c_,w(  tg  dB*  cosh  fcfo,,  +  z 

2  \    a  dx(       cosh  kh  a  dx{        cosh  Jfc/i 

'    b     lf-'g*,B,inh*(*'>  +  «)c-fa,t  |  tgfc^sinhfc^  +  z)  iut\ 

2  V     a            coshfc/i  a             cosh  kh       *     J 


XV 


jrp    _     9sinh2k(h0  +  z) 
<  4sinh2kh 


dB  nt      dB*  „ 


dx 


dx( 


(B.12) 
(B.13) 
(B.14) 


/: 


■dz  = 
9    d 


4dx( 
4dxe 


dxi  dx( 


dx 


dx. 


r'  sinh  2k(h0  +  z) 
JM         sinh  2kh 

I"  cosh 2k(h0  +  fj)  coah2k(h0  +  z) 

l2kBmh2k(h0  +  f})      2ksmh2k(h0  + fj). 


(B.15) 
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m  H  J_  \dBB.     gg  J  tn    (  co8h2k(h0  +  fj)  coah2k(h0  +  z)   \    , 

4dif[aif  5zf     J./-fcoV2*6inh2*(/i0  +  ij)      2kemh2k{h0  +  fj))  dz 

1  0  +  V)      4k*mnh2k{h0  +  fi)) 


pg  d 

4  dx( 

pg  d 
4  azf 

P9   d 
4  dx. 


3xf        aif   J  vj 


v2Asinh2Jk(A0  +  ^) 


/2khcoeh2kh         sinh  2fc/i    \ 
V4fc2sinh2Mi  "  4ifc2sinh2A:ftJ 

— 2  [2AAcoth2Jfc/i-  1] 


(B.16) 


Expanding  the  derivative  term  in  Eq.  (B.16) 


dxf 


dx{  dxi 


d2BB>+d-XB  +  2dBdB' 


dxi 


dx] 


dx:  dx{ 


(B.17) 


Recognizing  that  f  is  a  free  index  so  that  ■£-  =  Vy, 


d2B 
c 

d2B 


dxi 


=  V\B  =  -k2B 


dxi 


=  V2hB*  =  -k2B' 


which  results  from  use  of  the  Helmholtz  equation  which  governs  diffraction  in  water  of 
constant  depth  and 


dBdB* 

so  that  the  expanded  derivative  term  (B.17)  becomes 

32Bn,      d2B*  dBdB* 

^B   +  ~^TB  +  2irr  —  =  2pv\J3|'  -  2*2|B| 


(B.18) 


dxi 


dx2( 


dxi  dx( 


(B.19) 


and  the  result  is 


/-!  ('jf  ^^  dz=P{  (lV^!2  -  *2!fil)  W  [»*coth»A  - 


1]  (B.20) 


The  author  wishes  to  express  some  reservations  concerning  the  validity  of  Eq.  (B.20). 
The  shallow  water  asymptote  for  this  expression  is  zero  since  for  small  k(h  +  fj),  coth  2k{h  + 
*)  -  tk(k+fi)>  and  2k(h+fj)coth2k{h+fj)-l  S  0.  In  deep  water,  however,  coth2k(h+fj)  2  1 
so  that  2k(h  +  fj)  coth2k(h  +  ?)-lS  2kh.  If  the  derivative  terms  are  non-zero,  the  Szx 
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term  will  be  proportional  to  the  depth.  This  would  mean  that  deep  water  non-shoaling  and 
non-dissipative  waves  in  varying  depths  would  produce  gradients  in  the  radiation  stresses. 
For  plane  waves  this  is  not  a  problem,  since  the  complex  amplitude  terms  are  easily  shown 
to  be  zero.  This  is  seen  by  substitution  of  B  =  aeikx  and  B*  =  oe-'**  where  a=\B\. 
The  fourth  term  is  straight  forward 

«/  =  I  (- i£kB^hk(h0  +  z)     iut      igk     .sinhfc^  +  z)   ■  A 

2  V      o  cosh**       *        +    a  B        coshkh      e     )  <B-21) 

,    hj       1  g2k2  sinh2  k(h0  +  z)  /         .  ..  ,  ... 

(" }    =  i^       cosh*  I*        (2BB   -  BC         -  B*C      )  (R22) 

T^yy-  ^'^h2  *(*  +  *)„„„.     gksinh*k(h0  +  z)„^ 

{W  >        4  ^coshH(/l  +  ,-)2BB sinh2*ft         BB  (B-23) 

-  f\  Pw?dz  =  -P9k\B\>  r  **£%&££*  =  -eilB{>  U  _  ^h\  (B  U) 

'-*.  J -h0       smh2kh  41    '    \       sinh2JbA/     K        ' 

The  fifth  term  needs  no  evaluation  and  is  seen  to  cancel  exactly  the  second  term. 
The  sixth  and  last  term  of  Eq.  (B.6)  is  ^(^  and  is  evaluated  easily 

„'=i(Be—  +  BV-)  (B.25) 

(J)2  =  \  (Be"2-'  +  2BB*  +  B*e2"")  (B.26) 

W=^B*  =  M!  (B.27) 

fW=f|B|2  (B.28) 

which  cancels  part  of  the  result  from  the  fourth  term. 

Adding  the  above  results  for  all  the  terms  of  Eq.  (B.6)  yields 

c     _  ^  j^fdBdB^    1   /  2kh     \ 

5-'-T(*  (^^jJ^(1  +  sinh2^J  + 

6ap  H'd^bfc  +  dV^|2  "  *2^|2)2^(2^coth2^  -  1)]  }  (B.29) 


APPENDIX  C 
LAGRANGIAN  DERIVATION  OF  THE  GOVERNING  WAVE  EQUATION 


The  governing  equation  for  wave  motion  can  be  derived  through  use  of  a  variational 
principle  where  a  Lagragian  is  the  function  that  is  varied.  Luke  (1967)  showed  that  a 
Lagrangian  for  waves  is  given  by 

L = I-ko  r + ^(v*)2 + gz\ dz  (C1) 

Kirby  (1984)  used  a  time  averaged  Lagrangian  to  derive  a  wave  equation.  For  the  following 
derivation,  equation  (C.l)  will  be  used.  The  third  term  in  Eq.  (C.l)  is  determined  straight 
forward,  while  the  first  two  terms  need  to  be  expanded  about  fj  the  mean  water  line  in  a 
Taylor  series. 


2  2        J- 


4h  +  ^2  d2  +  („  _  n)4>t  lfi  +{n  _  n) Qj?  u        (c  2) 


ghl  ,  /"» 

ho 


2 
We  assume  that  r)  and  <f>  can  be  expressed  as 

i?    =    fj  +  erii 

<f>    =    <t>o  +  cf<f>i  (C.3) 

where  c  is  an  ordering  parameter  and  /  is  the  depth  dependency  as  defined  in  Eq.  (3.3).  For 
the  sake  of  simplicity  we  will  assume  that  the  gradient  of  fo  represents  the  steady  horizontal 
currents.  This  means  that  <f>0  is  not  a  function  of  time  and  that  W  will  not  be  included  as 
it  was  in  the  derivation  in  Chapter  3.  This  is  because  the  W  term  is  not  expanded  about 
*  =  fj  as  it  was  in  the  Green's  function  derivation  and  is  not  needed  in  order  to  obtain  the 
(V*  *  u)-£  term  which  as  will  be  seen  is  obtained  more  directly  in  the  present  derivation. 
Including  W  here  will  only  produce  terms  that  are  of  higher  order  in  the  mean  current. 
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Introducing  (C.3)  into  (C.2)  yields 


h^hhf  [fVkf  +  ffM]d*)  + 
j — h0  ) 

*1iK  +  «?,  |(V^o)    +  tVh<f>QVhfa  +  eVfc^x/,  |„  +0(es)l 
The  Lagrangian  1  is  assumed  to  be  of  the  form 


(C.4) 


L  =  Lo  +  eli  +  e2^2  +  eSI/3  H (C.5) 

It  then  follows  that 

r  _w2    gh20  ,  (v^0)2,L 

°~~2 2  2       (  °  +  ^  (C6) 

*a  =  9Hm  +  -j-*i,  +  y V^0 Vkh  [J"   Vhf  dz  +  J"    f,dz  \  +  m  (V^o)2        (C.7) 


and 


4>i?hfa  J_^  \nhf  +  f  f„]  dz  +  nifat  +  k v^v^j  +  mVh<f>ofaf,  |,  (c.8) 


To  obtain  the  desired  linear  governing  equation  L2  is  varied  with  respect  to  fa  and  then  L2 
is  varied  with  respect  to  r^  in  order  to  eliminate  ifc.  The  variational  principle  holds  that  the 
Lagrangian  is  constant  under  small  variations  of  the  unknown  parameters.  In  mathematical 
terms 

DL      « 

Dfa  =  °  <C'9) 

and 

Since  the  ordering  parameter  c  is  arbitrary  the  above  holds  for  all  L,-   I  =  0, 1, 2,  •  •  •.  In  the 
L2  term  fa  is  expressed  in  terms  of  its  derivatives.   This  requires  integration  by  parts  to 
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obtain  the  variation  with  respect  to  fa.  Upon  carrying  this  out  it  is  found  that 

Dfa    a<h    dt(d(<t>uy   Vfcla(v^1))  (C11) 

Applying  (C.ll)  to  (C.8)  and  recognizing  that  the  variation  dfa  is  arbitrary  yields  after 
carrying  out  the  integrations  in  L2 

-  Vfc(^Vkfc)  +  +i(°*~*CC')  -  nu  -  v*(vfcMi)  =  0  (C.12) 

which  is  rewritten  as 

9~Dt  +  W1  V**>  +  Vh(CCiVh<f>i)  ~  M°2  ~  k2CCg)  =  0  (C.13) 

Varying  L2  with  respect  to  rji  yields 

9m  +  <t>i,  +  VhfoVkti  =  0  (C.14) 

which  can  be  written  as 

lDfa 

Vl  =  —g-Dr  (Cl5) 

which  upon  substitution  into  (C.13)  yields  the  governing  equation 

D26i  ,       £><*, 

~5F  +  {Vh4>o)Dt   "  V^CC'VM  +  (**  ~  k2CCt)fa  =  0  (C.16) 

which  upon  making  the  substitution  that  if  =  Vfc#0  becomes 

«_  ft    . 

~W-  +  {Vh'a)^f-  V*(CC«V^i)  +  (**  -  *2<?C,)^  =  0  (C.17) 


APPENDIX  D 
OBTAINING  EQUATION  (3.79) 


In  this  Appendix  we  show  the  algebra  involved  in  obtaining  equation  (3.79)  from  equa- 
tion (3.51)  using  the  assumption  that  the  dispersion  relationship  is  written  as 

u  =  <r  +  kcosOU  m  i) 

We  start  by  using  the  following  form  of  $ 

JU -frdfelljUf  *-#*!  (D2) 

and  then  make  a  phase  shift  by  using  the  substitution 

A  =  A'ei(flco,!d'-fkco't<1*)  (D  3) 

so  that  in  the  final  form 

4,  =  -ig^ylei(fi^dz)  (D4) 

and  all  y  dependency  is  solely  in  the  complex  amplitude  A'.  For  reference  equation  (3.51) 
is  reproduced  below. 

(p4>z)x  +  {p4>v)v  +  k2p4>  +  [w2  -  a2  +  ,'w(Vfc  •£?)]£  +  iVw^ 
-{U24>z)x  -  {V2]>v)v  -  {UV4>X)V  -  (UVfa.  +  2iuU  •  Vh4>    =    0  (D.5) 

Writing  out  all  the  components  of  equation  (D.5)  and  rearranging  terms  gives 
(P  -  U2)t$x  +  (p  -  U2)^x  +  [{p  -  V2)4>y]v  +  k2pi  +  (u2  -  a2)l 
^Uai+wVw++iow$-(UVfa9-(UVfam+2iuU&  +  ^V49    =    0    (D.6) 
The  derivatives  of  <j>  defined  by  equation  D.2  are  as  follows 

&     =     -ig[(jj    +fJfccos^]c''(/tco"d*) 

120 


121 


4>xx 


-ig  [ (j\     +  2ik cosOj  +  « (jfc cos 6)x -  -  k2 cos2 $ -  j  e'Xf  *  eo"  <**) 
ik)y 


The  above  forms  of  the  derivatives  recognize  that  the  phase  function  e*(f kcottdx)  nas  y 
derivatives  since  *  =  k(x,y).  However  the  y  derivatives  of  the  phase  function  integral  are 
not  evaluated  since  the  phase  shift  (equation  D.3)  will  eliminate  any  y  dependence  from 
the  phase  function  and  place  all  y  dependance  in  the  complex  amplitude.    The  M) 

V  °  I XX 

term  is  neglected  following  the  assumption  explained  in  section  3.2.  Substituting  the  above 
derivatives  into  equation  (D.6)  gives 


(P  -  U2)x  (-)    e''(/4co''^)  +  ik  cos  6[p  -  U2)x-e^fkco'td') 
+2ikcoe6{p-U2)  (-)    eilfk™9d*)+i(kcoB0)x{p-U2)-eilfkco»dx) 

\& /  x  O 

-Jfc2  cos2  9  (p  -  U2)  -e'(/  *  co*  '  *)  +    (p-  V2]  (-c'(/  *  c°«  •  dzA  1 

+Jfe2p-e'(/*C0«'«'*)  +  (u2  -  a2)-ei(fkc'"'d*)  +  iu,Ut-ei<Jkco"d') 
o  'a  a 

+ivVv  -e'tf  * "■  '  d')  +  iwAe'lf  k cot  *  dx)  -  (UV)y  (-)    Jtf  k co» '  *0 

-  2UV  f  (£\    e'tf  *  co' ' dz)     -ik  cos  OUyV-e'lf  k "■ ' *) 

-t*costfJ7Vv-e,'(/*c»«'^)  -  i(*cos0)ytfV-e,'(/*co"<il) 

-2ikcoB0UV  (^eilfkc<*tdxA    -i[kca*8)vW±&! *"»•**) 

-{UV)x(^fkcottdzA    +2iu>u(-\    e«f*co,0dx) 
-2wkcoe9U-eif>fkcc*tdz)  +  2iuV  {-?(!*<*»'**)}      m    0(D7) 
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Dropping  terms  that  contain  products  of  derivatives  of  £  and  derivatives  of  the  current  as 
well  as  products  of  the  current  and  second  derivatives  of  £  (i.e.  terms  that  are  designated 
by  an  overbrace)  equation  (D.7)  is  rearranged  as 

{•  [k  cos  9 (p  -  V2)  +  uU]  x  j  +  2t  [k  cos  9 (p  -  U2)  +  uU]  (-)    +  iwA 

+  [-k2  cos2  9{p  -  U2)  +  k2p  +  u2-a2-2ujk  cos  01/]  -  -  i(k  cos  9)VUV- 

+i[Vv(cj  -  k  cos  9U)  -  {k  cos  9U)VV  -  k  cos  9UV]  -\  «<(/*««•««) 

+  K^*00"^)  1    +2ilV("-kco*eU)}(jei{!kC°'edx))       =    C(D.8) 
It  is  easily  shown  that 

kcos9(p-U2)  +  uU  =  a(Cgcos9  +  U)  (D.9) 

By  writing  the  dispersion  relationship  as 

a  =  u  -  kcos9U  (D.10) 

it  is  also  easily  shown  that 

u>2  -a2  =  2ukcos9U  -(kcos9U)2  (D.ll) 

From  equation  D.2  it  is  apparent  that  for  a  constant  ui, 

(kcos9)v  =  -Oy  (D.12) 

Thus  equation  (D.8),  becomes 

{•  \a{Cg  cos  9  +  U)]x  ^  +  2«r(C,  cos  9  +  U)  (-\ 
+  [k2p{l  -  cos2  9)]  £  +  iwA  +  tV,A  +  i*yV-\  jll *«••<*) 

+2,Va(^(/*-«^    +L^'(/»-'^l       =    o  (D.13) 


v      I    *-  'yjv 


where  a  i(fccos0),,C/Vj  term  has  been  dropped. 
At  this  point  a  phase  shift  is  introduced 


yj  _   At -«(/  k~  cot  f  dz-f  k cot  tdz) 


(D.14) 
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where  i  and  9  are  a  reference  wave  number  and  wave  angle.  For  computational  purposes  ifc 
will  be  an  averaged  wave  number  across  a  grid  row  and  9  will  be  taken  as  the  wave  angle  in 
the  absence  of  currents  using  Snell's  law.  Further  it  will  be  assumed  solely  for  computational 
purposes  that  $  is  equal  to  0.  This  gives 

i  [a(Cg cos 9  +  U)]x  ^  +  2ia(Cg cob9  +  U)  f-\ 
-2(Cg  cos  9  +  U)  cos  9{k-  k)A 

+N> -«•*')]  7  4®  J, 

+iwA  +  tV,A  +  2iVa(j)    +iavV-    =    0  (D.15) 

Expanding  the  derivatives  of  £  and  with  a  few  minor  manipulations  of  the  terms  the 
governing  equation  (3.79)  is  obtained. 


A 

X 


[C,  cos 9  +  U)AX  +  %{C,  cos 9  +  U) cos 9{k  -  k)A  +  -  \Ct*»9  +  U' 

2  [  a 

+VA"  +  I  (If)/  "  ^  -  «*'  "A  -  5  [CC-  (?)  J    +  f»    -    «D- "•) 
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particular  the  comparison  to  the  laboratory  experiments  of  Gourlay)  confirm  the  predictive 
capability  of  the  numerical  model.  The  simplifications  mentioned  above  make  this  possi- 
ble, for  there  must  be  a  mathematical  statement  of  the  physical  process  before  numerical 
modelling  of  the  process  can  even  be  attempted. 

In  making  recommendations  for  further  study  we  should  distinguish  between  short  range 
objectives  and  the  possibilities  that  may  become  open  to  numerical  modelers  by  advances 
in  the  computer  hardware  field.  It  is  not  difficult  to  envision  within  a  few  decades  or  less, 
super  computers  with  thousands  of  parallel  processors.  Such  advances  will  make  many  of  the 
present  computational  methods  obsolete.  For  example  the  parabolic  approximation  is  made 
solely  to  put  the  wave  equation  in  a  solvable  form  with  respect  to  the  capabilities  of  present 
day  computers.  With  sufficient  parallel  processors  it  is  conceivable  that  it  may  be  practical 
to  solve  the  full  elliptic  equation  explicitly.  We  shall  not  dwell  further  on  speculation  about 
future  generations  of  computers.  In  the  near  term,  the  numerical  model  presented  in  this 
report  can  be  refined  and  the  applications  can  be  extended.  The  model  is  constructed  in 
modular  form  so  that  any  new  formulation  for  any  of  the  component  processes  can  easily 
be  incorporated  into  the  model.  The  model  is  thus  open  for  continual  improvement. 

6.1     Improved  Wave  Formulations 

The  model  can  be  improved  by  better  wave  modeling  and  by  a  more  correct  formulation 
of  the  radiation  stress  terms  once  breaking  has  commenced.  The  formulation  used  in  this 
model  and  previous  models  for  the  radiation  stresses  initiates  set-up  and  longshore  currents 
at  the  breaker  line  whereas  as  explained  in  section  5.1  experiments  show  otherwise. 

An  improved  wave  model  will  assign  a  non-zero  value  of  the  amplitude  in  the  lee  of 
a  shore  parallel  breakwater.  Research  should  be  done  to  establish  a  method  of  doing  this 
so  as  to  enable  more  wave  energy  to  propagate  into  the  lee  of  the  structure.  Dalrymple 
and  Kirby  (1988)  have  found  a  method  to  determine  a  non-zero  wave  height  in  the  lee  of  a 
breakwater.  However  their  solution  is  in  the  form  of  a  Fourier  integral  which  is  not  suited 
to  the  complex  amplitude  form  required  in  the  present  model. 
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6.2     Including  Wave  Reflection 

The  author  has  made  unsuccessful  attempts  to  include  wave  reflection  in  the  model. 
Kirby  (1986)  used  coupled  equations  governing  the  interaction  of  a  forward  and  a  backward 
(reflected)  propagating  wave  that  are  coupled  through  a  slope  interaction.  He  was  able  to 
obtain  the  reflection  and  transmission  coefficients  over  submerged  breakwaters  that  agree 
well  with  the  experiments  of  Seelig  (1980).  Once  the  arrays  of  the  complex  amplitudes  of 
the  two  waves  is  determined,  the  radiation  stress  terms  can  be  determined  from  equations 
that  have  been  derived  following  the  same  procedure  as  used  by  Mei  (1972)  but  defining  4> 
as  being  comprised  of  two  components.  In  fact  this  procedure  can  also  be  used  to  obtain 
the  radiation  stress  terms  for  any  number  of  components.  The  only  limit  being  the  amount 
of  computation  required  as  the  number  of  components  increases. 

Smith  (1987)  following  the  work  of  Miles  (1967)  for  wave  reflection  and  transmission 
over  a  step  derived  equations  for  the  complex  reflection  and  transmission  from  a  submerged 
plane  barrier.  An  attempt  was  made  using  this  formulation  to  obtain  the  circulation  and 
wave  pattern  from  a  submerged  breakwater  represented  as  a  plane  impermeable  barrier. 
However  problems  of  numerical  stability  were  encountered.  It  is  suspected  that  the  sharp 
gradients  between  the  region  with  the  presence  of  the  reflected  wave  and  regions  of  its 
absence  adjacent  to  the  barrier  produced  unrealistic  gradients  in  the  radiation  stress  terms. 
An  attempt  was  made  to  include  the  reflected  wave  for  the  case  of  an  emerged  plane  barrier 
with  the  same  unstable  result.  This  is  a  topic  for  further  investigation. 

63     Extending  the  Model  to  Include  Sediment  Transport 

The  model  as  it  currently  stands  can  be  extended  to  include  sediment  transport.  The 
wave  model  can  be  slightly  amended  to  retain  wave  energy  dissipation  as  an  array.  If  a 
formulation  of  sediment  entrapment  can  be  made  in  terms  of  the  wave  energy  dissipation 
and  the  bottom  shear  stress,  everything  else  is  in  place  to  accomplish  the  sediment  transport 
model.  A  simple  numerical  mechanism  would  be  to  assume  steady  state  conditions  for 
certain  time  intervals  and  solve  for  the  wave  and  current  field  and  the  sediment  concentration 
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as  a  function  of  known  variables  and  then  solve  a  separate  model  for  sediment  balance  and 
then  update  the  depths  and  proceed  iterating  between  the  sediment  balance  model  and  the 
present  wave-current  model  marching  in  time.  Since  it  has  already  been  demonstrated  that 
the  wave-  current  model  gives  adequate  results  in  the  lee  of  a  shore  parallel  breakwater,  a 
good  test  for  such  a  sediment  transport  model  wold  be  to  see  if  it  could  yield  a  tombolo 
formation  in  the  lee  of  a  shore  parallel  breakwater. 


APPENDIX  A 
DERIVATION  OF  THE  DEPTH  AVERAGED  EQUATIONS 


In  this  appendix  a  derivation  of  the  governing  equations  for  the  circulation  model  is 
presented.  This  appendix  closely  follows  the  work  of  Ebersole  and  Dalrymple  (1979). 

Before  proceeding  with  the  derivations  of  the  depth  integrated  time  averaged  equations 
of  motion  and  the  continuity  equation  the  Leibnitz  Rule  for  taking  the  derivative  of  an 
integral  is  stated  since  it  is  invoked  throughout  the  derivation. 

The  Leibnitz  Rule  can  also  be  used  in  reverse  to  take  derivatives  outside  of  an  integral. 
The  kinematic  boundary  conditions  are  also  reproduced  since  the  use  of  the  Leibnitz  Rule 
will  yield  bottom  and  surface  terms  that  will  be  eliminated  through  use  of  these  boundary 
conditions.  On  the  free  surface  the  boundary  condition  is 

dr\  dr\  drj 

and  at  the  bottom  the  boundary  condition 

dh0  dh0  dh0 

In  the  circulation  model  it  is  assumed  that  the  still  water  depth  does  not  vary  with  time  so 
the  first  term  of  Eq.  (A.3)  is  dropped  (actually  for  time  varying  depth  essentially  the  same 
equations  are  obtained). 

The  continuity  equation  (or  the  conservation  of  mass  equation  as  it  is  also  called)  in 
three  dimensions  is 

gg  |   gtgO   ,   d(pv)      djpw) 


107 


108 
Integrating  over  the  depth  from  z  =  -k0  to  z  =  q  and  using  the  Leibnitz  Rule  to  take 
derivatives  outside  the  integral,  the  continuity  equation  becomes 


,     x   drj  dh0      d    f> 


fin  Ah  ^ 

-(H^-(H-A<,^  +  (^)n-(£uO_^    =    0  (A.5) 

The  terms  with  an  overbrace  are  p  times  the  free  surface  kinematic  boundary  condition 
and  the  terms  with  an  underbrace  are  p  times  the  bottom  boundary  condition.  Eliminating 
both  of  these  from  Eq.  (A.5)  gives 

Ft  U  pdz  +  rx  U  pudz+rv  U pvd2  =  0  (A6) 

Equation  (A.6)  is  now  time  averaged.    Time  averaging  means  averaging  over  a  wave 
period.  The  time  average  of  a  term  is  denoted  by  an  overbar.  Thus 

J=fIoFdt  <A-7) 

where  F  represents  any  term  and  T  is  the  wave  period.  Letting  u  and  v  be  comprised  of  a 
mean  flow,  a  wave  induced  flow,  and  an  arbitrary  or  random  turbulence  component 

u  =  u  +  ti  +  u'  (A.8) 

v  =  v  +  v  +  v'  (A.9) 

The  turbulent  fluctuations  are  assumed  to  be  of  a  high  frequency  relative  to  the  wave 
frequency  so  that  their  time  averages  over  the  wave  period  are  zero.  By  substituting  the 
above  expressions  for  u  and  v  into  Eq.  (A.6)  and  taking  time  averages  the  continuity  equation 
is 

Jl  {P(h0  +  fj)}  +  —  {pu(h0  +  fj)}  +  —  j*    ptdz+  (A.10) 


J-y  {*>(*.  +  *»  +  d-yf-hopidz  =  °  (A11) 
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The  time  averages  of  the  vertically  integrated  wave  induced  velocities  u  and  v  are  not  zero. 
These  quantities  are  known  as  the  wave  induced  mass  flux. 
Defining 

U  =  fi  +  u  and  V  =  v  +  v  (A. 12) 

where 


u  = 


(A.13) 


the  time  averaged  depth  integrated  continuity  equation  is 

Tt  +kU(K  +  f)  +  ^(*.  +  *)-0  (A.14) 

The  assumption  that  the  density  is  constant  in  both  time  and  space  has  used  in  obtaining 
both  Eq.  (A. 11)  and  Eq.  (A.14).  The  assumption  of  a  constant  bottom  has  also  been  used 
in  obtaining  (A.14).  If  the  bottom  is  not  constant  in  time  a  ^f  term  would  be  included  in 
Eq.  (A.14). 

The  x  and  y  momentum  equations  are  handled  in  the  same  manner  as  the  continuity 
equation  above.  The  x  and  y  direction  momentum  equations  are 

du      dtS       duv      duw  _  1  dp       1    drt *       3rvx       drMt , 

dt^  dx  +  dv  +  dz  ~  0di  +  ~o{~dx-  +  -a^  +  -fi7i  (A15) 


^      ^      M       dw       Idp      Idr^      dry        dr, 

dt*   dx   *  dy  +    dz         pdy  +  p{  dx    +1y~+~dT}  (A16) 

The  same  procedure  is  followed  for  both  the  x  and  y  equations.  Only  Eq.  (A.  15)  will  be 
dealt  with  here.  Integrating  the  left  hand  side  over  the  depth  and  applying  the  Leibnitz 
Rule  to  take  the  derivative  outside  the  integral  gives 

(LHS)   "   dll-Cdz~u«Tt-uj^jl+ 

J -h0  ''dx        -*•  dx 

(uw),,  -  {uw).ho  (A.17) 


a_ 

dx 
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The  terms  in  Eq.  (A.17)  that  are  designated  by  an  overbrace  are  u„  times  the  kinematic 
free  surface  condition  and  those  designated  by  an  underbrace  are  u_fco  times  the  bottom 
boundary  condition.  These  are  both  dropped  leaving 

(LHS)  =  §iJ"ho"**+£j"h/d*+§-J"hUvdz  (A.18) 

Substituting  the  definitions  given  for  u  and  v  in  Eqn.  (A.8)  and  (A.9)  and  recognizing  that 
the  turbulent  fluctuating  components  are  of  a  high  frequency  relative  to  the  wave  frequency 
so  that  the  time  averages  of  integrals  involving  the  product  of  a  turbulent  component  and 
a  wave  induced  component  are  zero  gives 

LHs "  HZ. s  * + mil.  •  * + j-XK  •** + iLLf"*    <ai9> 

+  Yxihy^dz  +  aiU  2™dz  +  aiU  ™dz  +  aiU  ™dz         (A20) 

B  ~H>    "  d  ~pi  d  ~n 

+  TyL.  U€  ^  +  TvU  » dz  +  TyU  «'*' dz  (A21) 

using  the  definitions  of  Eqs.  (A.12)  and  (A.13)  this  becomes 


LHS    =    |^  +  ,,  +  A^  +  ,-)+^^_^[_i_{/_:/i<(2} 


+hm^ « +fvi:. « *  -  £  y^ij^z? * 


(A.22) 


where  the  two  terms  designated  by  underbraces  result  from  being  added  to  complete  the 
U2  and  UV  terms. 

The  right  hand  side  after  integrating  over  the  depth  assuming  a  constant  density  and 
assuming  an  inviscid  fluid  such  that  no  horizontal  viscous  stresses  exist,  i.e.  r«  and  ryx  are 
zero,  is 

p  J-ho  dx      T  p  J.Ke   dz  aZ  ^2Z) 

Using  the  Leibnitz  Rule  and  averaging  over  a  wave  period  the  right  hand  side  becomes 


p  dx].ho  P      T  pP"dx  +  pP~h°  dx  +  pT""       pT"-»c 


(A.24) 


Ill 

r^"  and  r„_ko  are  the  time  averaged  surface  and  bottom  shear  stresses.  Assuming  that 
the  free  surface  pressure  p„  is  zero  (the  dynamic  free  surface  boundary  condition)  the  right 
hand  side  is 


RHS"  -pTxUpdz  +  -/-»•  *T+  ~^~  -^  (A.25) 

The  mean  pressure  at  the  bottom  can  be  expressed  as  the  sum  of  the  wave  induced  or 
dynamic  pressure  and  the  hydrostatic  pressure 

P^ho  =  Pdyn.ko  +  pg{h0  +  fj)  (A.26) 

so  that  the  product  p^^fe  can  be  manipulated  to  read 

v=s^k  -  \^^d~t + \ih{9{ho + n),}  ~  ^ +^Yx       (A-27> 

Terms  from  both  the  right  and  left  hand  sides  are  combined  as  the  radiation  stresses. 
Defining 


s" ~ I -J* + p"2) dz " 2pg(h° + ^ ~ pjir+^itf-  p<ldz^      (A-28) 

S,v  =  f\utdz  -  -^J^J^Z  (A.29) 

Sw  =  J_hy+f°2) dz  -  2P9^ho  +  *>?  ~  nr+  -) ij    Pidz)2       (A-3°) 

and  making  the  assumptions  that  the  product  of  the  dynamic  pressure  at  the  bottom  and 
the  bottom  slope  is  negligible,  that  the  &f2ko(u')*  dz  term  can  be  ignored  and  that  the 
lateral  friction  rt  caused  by  momentum  fluxes  due  to  turbulent  fluctuations  is  independent 
of  depth  and  that  r,  is  defined  by  r,  =  -pu'v',  the  x  direction  equation  of  momentum  is 
written  as 

*U(ha  +  If)  +  §-U\h0  +  fj)  +  ±UV(h0  +  fj)  +  l*p*  +  l*fJ3L 
at  ox  dy  p   dx        p    dy 

,  _.dfj      {h0  +  fj)dTt       1  1 

+,(*.  +  „)_  +  i-^-i^  -  -r„  +  -*.  =  0  (A.31) 

By  subtracting  the  depth  integrated  continuity  equation  and  dividing  by  the  total  depth, 
the  final  form  of  the  x  momentum  equation  is 

^  +  u^  +  vJy-  +  pJ^TT)[-aT  +  -af}  (A32) 
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+gdx  +  pd^-  7U^TT)T"  +  7{KTW)n*  =  °  (A-33) 

The  same  procedure  gives  as  the  depth  integrated  time  averaged  y  momentum  equation 

^  +  u—  +  VdV  +         ?         \dS">  ^  dS*v1 
dt+UdX+V^  +  rt^r)[-d?  +  -e?\  (A.34) 

+gdy +-p*;-  pv^T)r-« + pjh^nrf* = °  (A-35) 


APPENDIX  B 
DERIVATION  OF  THE  RADIATION  STRESS  TERMS 


The  basis  of  the  derivation  is  to  use  the  definition  of  the  velocity  potential  to  substitute 
for  the  wave  induced  velocities  in  the  definition  for  the  radiation  stresses. 

d<j> 


«i- 


dxp 


,       8<f> 

w=Tz 


where  (p  is  expressed  in  terms  of  the  general  complex  amplitude  function  B 


as 


where  ex.  represents  the  complex  conjugate  term.  Also  r}'  is  expressed  as 

rfml  (Be"*-*  +  c.c) 

where  the  prime  is  used  to  denote  wave  induced  quantities. 

For  convenience  Eq.  (2.50)  of  Section  2.2  is  reproduced  below. 

8+  =  pf^  u^dz  +  SaP  [J%dz  -P{(h0  +  fi?  +  f  W] 
An  expression  for  p  is  obtained  in  Section  2.2  and  is  reproduced  bel 


ow 


P[x,y,z)  =  pg(r)-z)  +  p        —L—dz-p(w>y 

Jm      uX( 

Substituting  Eq.  (B.5)  into  Eq.  (B.4)  gives 


Sap    =    P  u'au'0dz  +  6af) 


(B.l) 


(B.2) 


(B.3) 


(B.4) 


(B.5) 


-/_fc/K)2^-f(/»o  +  ^  +  fw 


(B.6) 
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The  terms  of  Eq.  (B.6)  are  now  evaluated  term  by  term. 
Term  1: 


a»  =  I  (_*9dB  cosh k(h0  +  z)      ■  ,      ig3B' cosh  *(*.  +  *)  M\ 
"      2\    aflia       coshjfcfc  ^  a  dxa        c^sh**       *     J  (BJ) 


i  mi  _  g2  cosh2  k(h0  +  z) 
4<r2  cosh2  kh 


«««/>  = 


dxa  dxp       dx*  dip      dxa  dxp  dxa  dxfi  * 


(B.8) 


pW    =     P9  cosh2  k{h0  +  z) 
a  fi  2Jfcsinh2A:/i 

pg  cosh2  k(h0  +  z) 
k  sinh  2kh 


dB  dB*      dB"  dB 


dxa  dip       dxa  dip 
dB  dB 


SR 


dxa  dip 


; 


(B.9) 


/      pu'au'pdz    = 

J  —h0 


ksinh2kh 

99 

ksmh.2kh 


'(  dB  dB*\\   ft  , 

=  nm(dBBB*\  i  /       abfc  \ 

4      \dxadx{,J  k2  V        sinh  2**/ 


(B.10) 


The  second  term: 


y_fc  M*? -*)<**=  pg(^-y) 


—  w — « — 


-Ac 


(B.ll) 


The  third  term: 


)eiu,t 


u'     B     1  /    «g  aB  cosh  fc(ft0  +  z)  c_,w(  tg  dB*  cosh  fcfo,,  +  z 

2  \    a  dx(       cosh  kh  a  dx{        cosh  Jfc/i 

'    b     lf-'g*,B,inh*(*'>  +  «)c-fa,t  |  tgfc^sinhfc^  +  z)  iut\ 

2  V     a            coshfc/i  a             cosh  kh       *     J 


XV 


jrp    _     9sinh2k(h0  +  z) 
<  4sinh2kh 


dB  nt      dB*  „ 


dx 


dx( 


(B.12) 
(B.13) 
(B.14) 


/: 


■dz  = 
9    d 


4dx( 
4dxe 


dxi  dx( 


dx 


dx. 


r'  sinh  2k(h0  +  z) 
JM         sinh  2kh 

I"  cosh 2k(h0  +  fj)  coah2k(h0  +  z) 

l2kBmh2k(h0  +  f})      2ksmh2k(h0  + fj). 


(B.15) 
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m  H  J_  \dBB.     gg  J  tn    (  co8h2k(h0  +  fj)  coah2k(h0  +  z)   \    , 

4dif[aif  5zf     J./-fcoV2*6inh2*(/i0  +  ij)      2kemh2k{h0  +  fj))  dz 

1  0  +  V)      4k*mnh2k{h0  +  fi)) 


pg  d 

4  dx( 

pg  d 
4  azf 

P9   d 
4  dx. 


3xf        aif   J  vj 


v2Asinh2Jk(A0  +  ^) 


/2khcoeh2kh         sinh  2fc/i    \ 
V4fc2sinh2Mi  "  4ifc2sinh2A:ftJ 

— 2  [2AAcoth2Jfc/i-  1] 


(B.16) 


Expanding  the  derivative  term  in  Eq.  (B.16) 


dxf 


dx{  dxi 


d2BB>+d-XB  +  2dBdB' 


dxi 


dx] 


dx:  dx{ 


(B.17) 


Recognizing  that  f  is  a  free  index  so  that  ■£-  =  Vy, 


d2B 
c 

d2B 


dxi 


=  V\B  =  -k2B 


dxi 


=  V2hB*  =  -k2B' 


which  results  from  use  of  the  Helmholtz  equation  which  governs  diffraction  in  water  of 
constant  depth  and 


dBdB* 

so  that  the  expanded  derivative  term  (B.17)  becomes 

32Bn,      d2B*  dBdB* 

^B   +  ~^TB  +  2irr  —  =  2pv\J3|'  -  2*2|B| 


(B.18) 


dxi 


dx2( 


dxi  dx( 


(B.19) 


and  the  result  is 


/-!  ('jf  ^^  dz=P{  (lV^!2  -  *2!fil)  W  [»*coth»A  - 


1]  (B.20) 


The  author  wishes  to  express  some  reservations  concerning  the  validity  of  Eq.  (B.20). 
The  shallow  water  asymptote  for  this  expression  is  zero  since  for  small  k(h  +  fj),  coth  2k{h  + 
*)  -  tk(k+fi)>  and  2k(h+fj)coth2k{h+fj)-l  S  0.  In  deep  water,  however,  coth2k(h+fj)  2  1 
so  that  2k(h  +  fj)  coth2k(h  +  ?)-lS  2kh.  If  the  derivative  terms  are  non-zero,  the  Szx 
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term  will  be  proportional  to  the  depth.  This  would  mean  that  deep  water  non-shoaling  and 
non-dissipative  waves  in  varying  depths  would  produce  gradients  in  the  radiation  stresses. 
For  plane  waves  this  is  not  a  problem,  since  the  complex  amplitude  terms  are  easily  shown 
to  be  zero.  This  is  seen  by  substitution  of  B  =  aeikx  and  B*  =  oe-'**  where  a=\B\. 
The  fourth  term  is  straight  forward 

«/  =  I  (- i£kB^hk(h0  +  z)     iut      igk     .sinhfc^  +  z)   ■  A 

2  V      o  cosh**       *        +    a  B        coshkh      e     )  <B-21) 

,    hj       1  g2k2  sinh2  k(h0  +  z)  /         .  ..  ,  ... 

(" }    =  i^       cosh*  I*        (2BB   -  BC         -  B*C      )  (R22) 

T^yy-  ^'^h2  *(*  +  *)„„„.     gksinh*k(h0  +  z)„^ 

{W  >        4  ^coshH(/l  +  ,-)2BB sinh2*ft         BB  (B-23) 

-  f\  Pw?dz  =  -P9k\B\>  r  **£%&££*  =  -eilB{>  U  _  ^h\  (B  U) 

'-*.  J -h0       smh2kh  41    '    \       sinh2JbA/     K        ' 

The  fifth  term  needs  no  evaluation  and  is  seen  to  cancel  exactly  the  second  term. 
The  sixth  and  last  term  of  Eq.  (B.6)  is  ^(^  and  is  evaluated  easily 

„'=i(Be—  +  BV-)  (B.25) 

(J)2  =  \  (Be"2-'  +  2BB*  +  B*e2"")  (B.26) 

W=^B*  =  M!  (B.27) 

fW=f|B|2  (B.28) 

which  cancels  part  of  the  result  from  the  fourth  term. 

Adding  the  above  results  for  all  the  terms  of  Eq.  (B.6)  yields 

c     _  ^  j^fdBdB^    1   /  2kh     \ 

5-'-T(*  (^^jJ^(1  +  sinh2^J  + 

6ap  H'd^bfc  +  dV^|2  "  *2^|2)2^(2^coth2^  -  1)]  }  (B.29) 


APPENDIX  C 
LAGRANGIAN  DERIVATION  OF  THE  GOVERNING  WAVE  EQUATION 


The  governing  equation  for  wave  motion  can  be  derived  through  use  of  a  variational 
principle  where  a  Lagragian  is  the  function  that  is  varied.  Luke  (1967)  showed  that  a 
Lagrangian  for  waves  is  given  by 

L = I-ko  r + ^(v*)2 + gz\ dz  (C1) 

Kirby  (1984)  used  a  time  averaged  Lagrangian  to  derive  a  wave  equation.  For  the  following 
derivation,  equation  (C.l)  will  be  used.  The  third  term  in  Eq.  (C.l)  is  determined  straight 
forward,  while  the  first  two  terms  need  to  be  expanded  about  fj  the  mean  water  line  in  a 
Taylor  series. 


2  2        J- 


4h  +  ^2  d2  +  („  _  n)4>t  lfi  +{n  _  n) Qj?  u        (c  2) 


ghl  ,  /"» 

ho 


2 
We  assume  that  r)  and  <f>  can  be  expressed  as 

i?    =    fj  +  erii 

<f>    =    <t>o  +  cf<f>i  (C.3) 

where  c  is  an  ordering  parameter  and  /  is  the  depth  dependency  as  defined  in  Eq.  (3.3).  For 
the  sake  of  simplicity  we  will  assume  that  the  gradient  of  fo  represents  the  steady  horizontal 
currents.  This  means  that  <f>0  is  not  a  function  of  time  and  that  W  will  not  be  included  as 
it  was  in  the  derivation  in  Chapter  3.  This  is  because  the  W  term  is  not  expanded  about 
*  =  fj  as  it  was  in  the  Green's  function  derivation  and  is  not  needed  in  order  to  obtain  the 
(V*  *  u)-£  term  which  as  will  be  seen  is  obtained  more  directly  in  the  present  derivation. 
Including  W  here  will  only  produce  terms  that  are  of  higher  order  in  the  mean  current. 
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Introducing  (C.3)  into  (C.2)  yields 


h^hhf  [fVkf  +  ffM]d*)  + 
j — h0  ) 

*1iK  +  «?,  |(V^o)    +  tVh<f>QVhfa  +  eVfc^x/,  |„  +0(es)l 
The  Lagrangian  1  is  assumed  to  be  of  the  form 


(C.4) 


L  =  Lo  +  eli  +  e2^2  +  eSI/3  H (C.5) 

It  then  follows  that 

r  _w2    gh20  ,  (v^0)2,L 

°~~2 2  2       (  °  +  ^  (C6) 

*a  =  9Hm  +  -j-*i,  +  y V^0 Vkh  [J"   Vhf  dz  +  J"    f,dz  \  +  m  (V^o)2        (C.7) 


and 


4>i?hfa  J_^  \nhf  +  f  f„]  dz  +  nifat  +  k v^v^j  +  mVh<f>ofaf,  |,  (c.8) 


To  obtain  the  desired  linear  governing  equation  L2  is  varied  with  respect  to  fa  and  then  L2 
is  varied  with  respect  to  r^  in  order  to  eliminate  ifc.  The  variational  principle  holds  that  the 
Lagrangian  is  constant  under  small  variations  of  the  unknown  parameters.  In  mathematical 
terms 

DL      « 

Dfa  =  °  <C'9) 

and 

Since  the  ordering  parameter  c  is  arbitrary  the  above  holds  for  all  L,-   I  =  0, 1, 2,  •  •  •.  In  the 
L2  term  fa  is  expressed  in  terms  of  its  derivatives.   This  requires  integration  by  parts  to 
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obtain  the  variation  with  respect  to  fa.  Upon  carrying  this  out  it  is  found  that 

Dfa    a<h    dt(d(<t>uy   Vfcla(v^1))  (C11) 

Applying  (C.ll)  to  (C.8)  and  recognizing  that  the  variation  dfa  is  arbitrary  yields  after 
carrying  out  the  integrations  in  L2 

-  Vfc(^Vkfc)  +  +i(°*~*CC')  -  nu  -  v*(vfcMi)  =  0  (C.12) 

which  is  rewritten  as 

9~Dt  +  W1  V**>  +  Vh(CCiVh<f>i)  ~  M°2  ~  k2CCg)  =  0  (C.13) 

Varying  L2  with  respect  to  rji  yields 

9m  +  <t>i,  +  VhfoVkti  =  0  (C.14) 

which  can  be  written  as 

lDfa 

Vl  =  —g-Dr  (Cl5) 

which  upon  substitution  into  (C.13)  yields  the  governing  equation 

D26i  ,       £><*, 

~5F  +  {Vh4>o)Dt   "  V^CC'VM  +  (**  ~  k2CCt)fa  =  0  (C.16) 

which  upon  making  the  substitution  that  if  =  Vfc#0  becomes 

«_  ft    . 

~W-  +  {Vh'a)^f-  V*(CC«V^i)  +  (**  -  *2<?C,)^  =  0  (C.17) 


APPENDIX  D 
OBTAINING  EQUATION  (3.79) 


In  this  Appendix  we  show  the  algebra  involved  in  obtaining  equation  (3.79)  from  equa- 
tion (3.51)  using  the  assumption  that  the  dispersion  relationship  is  written  as 

u  =  <r  +  kcosOU  m  i) 

We  start  by  using  the  following  form  of  $ 

JU -frdfelljUf  *-#*!  (D2) 

and  then  make  a  phase  shift  by  using  the  substitution 

A  =  A'ei(flco,!d'-fkco't<1*)  (D  3) 

so  that  in  the  final  form 

4,  =  -ig^ylei(fi^dz)  (D4) 

and  all  y  dependency  is  solely  in  the  complex  amplitude  A'.  For  reference  equation  (3.51) 
is  reproduced  below. 

(p4>z)x  +  {p4>v)v  +  k2p4>  +  [w2  -  a2  +  ,'w(Vfc  •£?)]£  +  iVw^ 
-{U24>z)x  -  {V2]>v)v  -  {UV4>X)V  -  (UVfa.  +  2iuU  •  Vh4>    =    0  (D.5) 

Writing  out  all  the  components  of  equation  (D.5)  and  rearranging  terms  gives 
(P  -  U2)t$x  +  (p  -  U2)^x  +  [{p  -  V2)4>y]v  +  k2pi  +  (u2  -  a2)l 
^Uai+wVw++iow$-(UVfa9-(UVfam+2iuU&  +  ^V49    =    0    (D.6) 
The  derivatives  of  <j>  defined  by  equation  D.2  are  as  follows 

&     =     -ig[(jj    +fJfccos^]c''(/tco"d*) 
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4>xx 


-ig  [ (j\     +  2ik cosOj  +  « (jfc cos 6)x -  -  k2 cos2 $ -  j  e'Xf  *  eo"  <**) 
ik)y 


The  above  forms  of  the  derivatives  recognize  that  the  phase  function  e*(f kcottdx)  nas  y 
derivatives  since  *  =  k(x,y).  However  the  y  derivatives  of  the  phase  function  integral  are 
not  evaluated  since  the  phase  shift  (equation  D.3)  will  eliminate  any  y  dependence  from 
the  phase  function  and  place  all  y  dependance  in  the  complex  amplitude.    The  M) 

V  °  I XX 

term  is  neglected  following  the  assumption  explained  in  section  3.2.  Substituting  the  above 
derivatives  into  equation  (D.6)  gives 


(P  -  U2)x  (-)    e''(/4co''^)  +  ik  cos  6[p  -  U2)x-e^fkco'td') 
+2ikcoe6{p-U2)  (-)    eilfk™9d*)+i(kcoB0)x{p-U2)-eilfkco»dx) 

\& /  x  O 

-Jfc2  cos2  9  (p  -  U2)  -e'(/  *  co*  '  *)  +    (p-  V2]  (-c'(/  *  c°«  •  dzA  1 

+Jfe2p-e'(/*C0«'«'*)  +  (u2  -  a2)-ei(fkc'"'d*)  +  iu,Ut-ei<Jkco"d') 
o  'a  a 

+ivVv  -e'tf  * "■  '  d')  +  iwAe'lf  k cot  *  dx)  -  (UV)y  (-)    Jtf  k co» '  *0 

-  2UV  f  (£\    e'tf  *  co' ' dz)     -ik  cos  OUyV-e'lf  k "■ ' *) 

-t*costfJ7Vv-e,'(/*c»«'^)  -  i(*cos0)ytfV-e,'(/*co"<il) 

-2ikcoB0UV  (^eilfkc<*tdxA    -i[kca*8)vW±&! *"»•**) 

-{UV)x(^fkcottdzA    +2iu>u(-\    e«f*co,0dx) 
-2wkcoe9U-eif>fkcc*tdz)  +  2iuV  {-?(!*<*»'**)}      m    0(D7) 
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Dropping  terms  that  contain  products  of  derivatives  of  £  and  derivatives  of  the  current  as 
well  as  products  of  the  current  and  second  derivatives  of  £  (i.e.  terms  that  are  designated 
by  an  overbrace)  equation  (D.7)  is  rearranged  as 

{•  [k  cos  9 (p  -  V2)  +  uU]  x  j  +  2t  [k  cos  9 (p  -  U2)  +  uU]  (-)    +  iwA 

+  [-k2  cos2  9{p  -  U2)  +  k2p  +  u2-a2-2ujk  cos  01/]  -  -  i(k  cos  9)VUV- 

+i[Vv(cj  -  k  cos  9U)  -  {k  cos  9U)VV  -  k  cos  9UV]  -\  «<(/*««•««) 

+  K^*00"^)  1    +2ilV("-kco*eU)}(jei{!kC°'edx))       =    C(D.8) 
It  is  easily  shown  that 

kcos9(p-U2)  +  uU  =  a(Cgcos9  +  U)  (D.9) 

By  writing  the  dispersion  relationship  as 

a  =  u  -  kcos9U  (D.10) 

it  is  also  easily  shown  that 

u>2  -a2  =  2ukcos9U  -(kcos9U)2  (D.ll) 

From  equation  D.2  it  is  apparent  that  for  a  constant  ui, 

(kcos9)v  =  -Oy  (D.12) 

Thus  equation  (D.8),  becomes 

{•  \a{Cg  cos  9  +  U)]x  ^  +  2«r(C,  cos  9  +  U)  (-\ 
+  [k2p{l  -  cos2  9)]  £  +  iwA  +  tV,A  +  i*yV-\  jll *«••<*) 

+2,Va(^(/*-«^    +L^'(/»-'^l       =    o  (D.13) 


v      I    *-  'yjv 


where  a  i(fccos0),,C/Vj  term  has  been  dropped. 
At  this  point  a  phase  shift  is  introduced 


yj  _   At -«(/  k~  cot  f  dz-f  k cot  tdz) 


(D.14) 
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where  i  and  9  are  a  reference  wave  number  and  wave  angle.  For  computational  purposes  ifc 
will  be  an  averaged  wave  number  across  a  grid  row  and  9  will  be  taken  as  the  wave  angle  in 
the  absence  of  currents  using  Snell's  law.  Further  it  will  be  assumed  solely  for  computational 
purposes  that  $  is  equal  to  0.  This  gives 

i  [a(Cg cos 9  +  U)]x  ^  +  2ia(Cg cob9  +  U)  f-\ 
-2(Cg  cos  9  +  U)  cos  9{k-  k)A 

+N> -«•*')]  7  4®  J, 

+iwA  +  tV,A  +  2iVa(j)    +iavV-    =    0  (D.15) 

Expanding  the  derivatives  of  £  and  with  a  few  minor  manipulations  of  the  terms  the 
governing  equation  (3.79)  is  obtained. 


A 

X 


[C,  cos 9  +  U)AX  +  %{C,  cos 9  +  U) cos 9{k  -  k)A  +  -  \Ct*»9  +  U' 

2  [  a 

+VA"  +  I  (If)/  "  ^  -  «*'  "A  -  5  [CC-  (?)  J    +  f»    -    «D- "•) 
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